export Quad

mutable struct Quad <: AbstractElement
    id::String
    hid::Int

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

    elm_type::String
    mass_type::String

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

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

    vertices=[n.loc for n in (node1,node2,node3,node4)]
    #check aspect ratio
    a=vertices[1]-vertices[3]
    b=vertices[2]-vertices[4]
    if !(0.1<norm(a)/norm(b)<10)
        @warn "Quad "*string(id)*" has an aspect ratio larger than 10, may lead to element illness."
    end
    #check wrong sequence
    for i in 1:4
        j=i%4+1
        if (vertices[i]-o)×(vertices[j]-o)⋅csys.T[:,3]<0
            @warn  "Quad "*string(id)*" has a wrong sequence of nodes, may lead to element illness."
            break
        end
    end
    #check devolution vertex
    for i in 1:4
        j=i%4+1
        if norm(vertices[i]-vertices[j])<1e-9
            @warn  "Quad "*string(id)*" has only 3 independent nodes, may lead to element illness."
            break
        end
    end
    #check distorsion
    for i in 1:4
        v=vertices[i]-o
        tol=cos(5*π/180)
        if abs(v⋅csys.T[:,3])/norm(csys.T[:,3])/norm(v)>tol
            @warn  "Quad "*string(id)*" is too distorted, may lead to element illness."
            break
        end
    end
    #check convexness
    for i in 1:4
        j=(i+1)%4+1
        k=i%4+1
        if norm((vertices[i]+vertices[j])/2-o)>=norm(vertices[k]-o)
            @warn  "Quad "*string(id)*" is not convex, may lead to element illness."
            break
        end
    end

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

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

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

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

function integrateK(elm::Quad)::SparseMatrixCSC{Float64}
    if elm.elm_type=="GQ12"
        K=K_GQ12(elm)
    elseif elm.elm_type=="TMQ"
        K=K_TMQ(elm)
    elseif elm.elm_type=="DKQ"
        K=K_DKQ(elm)
    elseif elm.elm_type=="DKGQ"
        K=K_DKGQ(elm)
    elseif elm.elm_type=="TMGQ"
        K=K_TMGQ(elm)
    else
        throw("Quad element type error!")
    end
    return K
end

function integrateKσ(elm::Quad,σ)::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)
    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


    function GtΣG(ξ,η)
        σ[1,1]= G*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2) + (1/2)*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2)*λ
        σ[2,2]= G*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2) + (1/2)*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2)*λ

        Σᵉ=[σ[1,1] 0;
            0 σ[2,2]]

        G=       [0  (-1/4)*(1 - η)              0  (-1/4)*(1 + η)               0  (1/4)*(1 - η)              0  (1/4)*(1 + η);
     (-1/4)*(1 - ξ)               0  (1/4)*(1 - ξ)               0  (-1/4)*(1 + ξ)              0  (1/4)*(1 + ξ)              0]

        J[1,1]=(-1/4)*(1 - η)*x₁ + (1/4)*(1 - η)*x₃ + (-1/4)*(1 + η)*x₂ + (1/4)*(1 + η)*x₄
        J[1,2]=(-1/4)*(1 - η)*y₂ + (1/4)*(1 - η)*y₃ + (-1/4)*(1 + η)*y₂ + (1/4)*(1 + η)*y₄
        J[2,1]=(-1/4)*(1 - ξ)*x₁ + (1/4)*(1 - ξ)*x₂ + (-1/4)*(1 + ξ)*x₃ + (1/4)*(1 + ξ)*x₄
        J[2,2]=(-1/4)*(1 + ξ)*y₃ + (1/4)*(1 + ξ)*y₄
        Gᵉ'*Σᵉ*Gᵉ*det(J)
    end
    Kσ=hquadrature(GtΣG,[-1,1],[1,1])
end

function integrateKu(elm::Quad,u::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)
    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 integrateK(elm,p_delta=false,σ=nothing,large_deform=false,u=nothing)
#     K=integrateK(elm)
#     if p_delta
#         K+=integrateKσ(elm,σ)
#     elseif large_deform
#         K+=integrateKu(elm,u)
#     end
#     return K
# end

function integrateM(elm::Quad)::SparseMatrixCSC{Float64}
    ρ=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)
    M=Matrix{Float64}(undef,24,24)
    J=Matrix{Float64}(undef,2,2)
    function ρNtN(x)
        ξ,η=x[1],x[2]
        # include(joinpath(dirname(@__FILE__),"jacobi.jl"))
        J[1,1]=(-1/4)*(1 - η)*x₁ + (1/4)*(1 - η)*x₃ + (-1/4)*(1 + η)*x₂ + (1/4)*(1 + η)*x₄
        J[1,2]=(-1/4)*(1 - η)*y₂ + (1/4)*(1 - η)*y₃ + (-1/4)*(1 + η)*y₂ + (1/4)*(1 + η)*y₄
        J[2,1]=(-1/4)*(1 - ξ)*x₁ + (1/4)*(1 - ξ)*x₂ + (-1/4)*(1 + ξ)*x₃ + (1/4)*(1 + ξ)*x₄
        J[2,2]=(-1/4)*(1 + ξ)*y₃ + (1/4)*(1 + ξ)*y₄

        M*det(J)
    end
    return hcubature(ρNtN,[-1,1],[1,1])[1]*t
end

function integrateP!(elm::Quad,elm_force)
    if elm.elm_type=="GQ12"
        P=P_GQ12(elm)
    elseif elm.elm_type=="TMQ"
        P=P_TMQ(elm)
    elseif elm.elm_type=="DKQ"
        P=P_DKQ(elm)
    elseif elm.elm_type=="DKGQ"
        P=P_DKGQ(elm)
    elseif elm.elm_type=="TMGQ"
        P=P_TMGQ(elm)
    else
        throw("Quad element type error!")
    end
    elm.Pᵉ=P
end
