#====================================================================
  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 findStrain!(mapDict::Dict{Int64, Int64}, Ïµ::Array{Float64, 1}, âˆ‚Ï•_âˆ‚x::Array{Float64,2},
    solAtNodes::Array{Float64, 1}, problemDim::Int64)

    fill!(Ïµ, 0.0)
    for a âˆˆ 1:size(âˆ‚Ï•_âˆ‚x, 1)
        for j âˆˆ 1:problemDim
            for i âˆˆ 1:j
                ij::Int64 = RapidFEM.getVoigtIndex(mapDict, i, j)
                c1::Float64 = (i==j) ? 0.5 : 1.0
                Ïµ[ij] += c1*(âˆ‚Ï•_âˆ‚x[a,j]*solAtNodes[problemDim*(a-1)+i]+âˆ‚Ï•_âˆ‚x[a,i]*solAtNodes[problemDim*(a-1)+j])
            end
        end
    end
end

function local_âˆ‡v_Cáµ€_âˆ‡u!(K::Array{Float64,2}, tensorMap_N_PlasticData::T,
    problemDim::Int64, element::AbstractElement, elementNo::Int64, shapeFunction::Array{ShapeFunction},
    coordArray::Array{Float64,2}; kwargs4function...)where T

    mapDict::Dict{Int64, Int64} = tensorMap_N_PlasticData[1]
    C::Array{Float64,2} = tensorMap_N_PlasticData[2]
    model::PlasticModel = tensorMap_N_PlasticData[3]
    modelParams::ModelParams  = tensorMap_N_PlasticData[4]
    stateDict  = tensorMap_N_PlasticData[5]
    stateDictBuffer  = tensorMap_N_PlasticData[6]
    #stateDictBufferCopy = deepcopy(stateDictBuffer)
    lastSoln::Array{Float64, 1} = tensorMap_N_PlasticData[7]
    plasticVars::PlasticVars = SmallStrainPlastic.initPlasticVars(model)
    plasticVars.C = C

    solAtNodes::Array{Float64, 1} = getSolAtElement(lastSoln, element, problemDim)
    âˆ‚Î¾_âˆ‚xFunc::Function = getFunction_âˆ‚Î¾_âˆ‚x(element)
    dÎ©Func::Function = getFunction_dÎ©(element)
    noOfIpPoints::Int64 = length(shapeFunction)
    noOfNodes::Int64 = size(shapeFunction[1].âˆ‚Ï•_âˆ‚Î¾,1)
    Ïµ::Array{Float64, 1} = zeros(model.ÏµVoigtSize)
    for ipNo::Int64 âˆˆ 1:noOfIpPoints
        âˆ‚x_âˆ‚Î¾::Array{Float64,2} = get_âˆ‚x_âˆ‚Î¾(coordArray, shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾)
        âˆ‚Î¾_dx::Array{Float64,2} = âˆ‚Î¾_âˆ‚xFunc(âˆ‚x_âˆ‚Î¾)
        x::Array{Float64, 1} = getInterpolated_x(coordArray, shapeFunction[ipNo].Ï•)
        dÎ©::Float64 = dÎ©Func(âˆ‚x_âˆ‚Î¾, shapeFunction[ipNo].ipData)
        âˆ‚Ï•_âˆ‚x::Array{Float64} = shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾*âˆ‚Î¾_dx
        findStrain!(mapDict, Ïµ, âˆ‚Ï•_âˆ‚x,  solAtNodes, problemDim)
        plasticVars.Ïµ = SmallStrainPlastic.getContinuumMandelStrain(Ïµ)

        #getState!(plasticVars.Ïµáµ–, plasticVars.Î±, stateDict, elementNo, ipNo)
        #plasticVars.Cáµ€ = SmallStrainPlastic.findNumerical_Cáµ€(plasticVars, model,
        #modelParams, stateDict,  elementNo, ipNo)
        #plasticVars.Cáµ€ = plasticVars.C
        SmallStrainPlastic.checkPlasticState!(plasticVars, model,
        modelParams, stateDict, stateDictBuffer,  elementNo, ipNo; algoTangent = true)
        #println("plasticVars.Ïµáµ– at âˆ‡v_Cáµ€_âˆ‡u= ", plasticVars.Ïµáµ–)
        #println("plasticVars.Cáµ€ = ", plasticVars.Cáµ€)
        for b::Int64 âˆˆ 1:noOfNodes
            for a::Int64 âˆˆ 1:noOfNodes
                for l::Int64 âˆˆ 1:problemDim
                    for k::Int64 âˆˆ 1:l
                        kl::Int64 = RapidFEM.getVoigtIndex(mapDict, k, l)
                        c2::Float64 = (k==l) ? 0.5 : 1.0
                        #c2 = 0.5
                        for j::Int64 âˆˆ 1:problemDim
                            for i::Int64 âˆˆ 1:j
                                ij::Int64 = RapidFEM.getVoigtIndex(mapDict, i, j)
                                c1::Float64 = (i==j) ? 0.5 : 1.0
                                #c1 = 0.5
                                K[problemDim*(a-1)+i,problemDim*(b-1)+k] += c1*c2*âˆ‚Ï•_âˆ‚x[a,j]*plasticVars.Cáµ€[ij,kl]*âˆ‚Ï•_âˆ‚x[b,l]*dÎ©
                                K[problemDim*(a-1)+j,problemDim*(b-1)+l] += c1*c2*âˆ‚Ï•_âˆ‚x[a,i]*plasticVars.Cáµ€[ij,kl]*âˆ‚Ï•_âˆ‚x[b,k]*dÎ©
                                K[problemDim*(a-1)+j,problemDim*(b-1)+k] += c1*c2*âˆ‚Ï•_âˆ‚x[a,i]*plasticVars.Cáµ€[ij,kl]*âˆ‚Ï•_âˆ‚x[b,l]*dÎ©
                                K[problemDim*(a-1)+i,problemDim*(b-1)+l] += c1*c2*âˆ‚Ï•_âˆ‚x[a,j]*plasticVars.Cáµ€[ij,kl]*âˆ‚Ï•_âˆ‚x[b,k]*dÎ©
                            end
                        end
                    end
                end
            end
        end
        plasticVars = SmallStrainPlastic.initPlasticVars(model)
        plasticVars.C = C
    end
    #SmallStrainPlastic.updateStateDict!(stateDictBuffer, stateDictBufferCopy)
    return nothing
end

function local_âˆ‡v_Ïƒ_Vector!(f::Vector, tensorMap_N_PlasticData::T, problemDim::Int64,
    element::AbstractElement, elementNo::Int64, shapeFunction::Array{ShapeFunction},
    coordArray::Array{Float64,2}; kwargs4function...)where T

    mapDict::Dict{Int64, Int64} = tensorMap_N_PlasticData[1]
    C::Array{Float64,2} = tensorMap_N_PlasticData[2]
    model::PlasticModel = tensorMap_N_PlasticData[3]
    modelParams::ModelParams  = tensorMap_N_PlasticData[4]
    stateDict  = tensorMap_N_PlasticData[5]
    stateDictBuffer  = tensorMap_N_PlasticData[6]
    #stateDictBufferCopy = deepcopy(stateDictBuffer)
    lastSoln::Array{Float64, 1} = tensorMap_N_PlasticData[7]
    plasticVars::PlasticVars = SmallStrainPlastic.initPlasticVars(model)
    plasticVars.C = C

    solAtNodes::Array{Float64, 1} = getSolAtElement(lastSoln, element, problemDim)
    âˆ‚Î¾_âˆ‚xFunc::Function = getFunction_âˆ‚Î¾_âˆ‚x(element)
    dÎ©Func::Function = getFunction_dÎ©(element)
    noOfIpPoints::Int64 = length(shapeFunction)
    noOfNodes::Int64 = size(shapeFunction[1].âˆ‚Ï•_âˆ‚Î¾,1)
    Ïµ::Array{Float64, 1} = zeros(model.ÏµVoigtSize)
    for ipNo âˆˆ 1:noOfIpPoints
        #Ïµ::Array{Float64, 1} = zeros(model.ÏµVoigtSize)
        âˆ‚x_âˆ‚Î¾::Array{Float64,2} = get_âˆ‚x_âˆ‚Î¾(coordArray, shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾)
        âˆ‚Î¾_dx::Array{Float64,2} = âˆ‚Î¾_âˆ‚xFunc(âˆ‚x_âˆ‚Î¾)
        dÎ©::Float64 = dÎ©Func(âˆ‚x_âˆ‚Î¾, shapeFunction[ipNo].ipData)
        Ï•::Array{Float64,1} = shapeFunction[ipNo].Ï•
        x::Array{Float64,1} = getInterpolated_x(coordArray, Ï•)
        âˆ‚Ï•_âˆ‚x::Array{Float64} = shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾*âˆ‚Î¾_dx
        findStrain!(mapDict, Ïµ, âˆ‚Ï•_âˆ‚x,  solAtNodes, problemDim)
        plasticVars.Ïµ = SmallStrainPlastic.getContinuumMandelStrain(Ïµ)
        #getState!(plasticVars.Ïµáµ–, plasticVars.Î±, stateDict, elementNo, ipNo)
        SmallStrainPlastic.checkPlasticState!(plasticVars, model,
        modelParams, stateDict, stateDictBuffer,  elementNo, ipNo)
        #if ipNo ==1
        #    println("plasticVars.Ïµ âˆ‡v_Ïƒ= \n", plasticVars.Ïµ)
        #    println("f = ", SmallStrainPlastic.ğ�’‡_j2(plasticVars.Ïƒ_voigt, plasticVars.q, plasticVars,modelParams), " plasticVars.Ïƒ_voigt =", plasticVars.Ïƒ_voigt)
        #end
        for a âˆˆ 1:noOfNodes
            for j::Int64 âˆˆ 1:problemDim
                for i::Int64 âˆˆ 1:j
                    ij::Int64 = RapidFEM.getVoigtIndex(mapDict, i, j)
                    c1::Float64 = (i==j) ? 0.5 : 1.0
                    f[problemDim*(a-1)+i] += c1*âˆ‚Ï•_âˆ‚x[a,j]*plasticVars.Ïƒ_voigt[ij]*dÎ©
                    f[problemDim*(a-1)+j] += c1*âˆ‚Ï•_âˆ‚x[a,i]*plasticVars.Ïƒ_voigt[ij]*dÎ©
                end
            end
        end
        plasticVars = SmallStrainPlastic.initPlasticVars(model)
        plasticVars.C = C
    end
    #SmallStrainPlastic.updateStateDict!(stateDictBuffer, stateDictBufferCopy)
    return nothing
end


function gaussian_Ïµáµ–(tensorMap_N_PlasticData::T,
    solAtNodes::Array{Float64,1}, problemDim::Int64,
    element::AbstractElement, elementNo::Int64, shapeFunction::Array{ShapeFunction},
    coordArray::Array{Float64,2}; kwargs4function...)::Array{Array{Float64,1},1} where T
    mapDict::Dict{Int64, Int64} = tensorMap_N_PlasticData[1]
    C::Array{Float64,2} = tensorMap_N_PlasticData[2]
    model::PlasticModel = tensorMap_N_PlasticData[3]
    modelParams::ModelParams  = tensorMap_N_PlasticData[4]
    stateDict  = tensorMap_N_PlasticData[5]
    stateDictBuffer  = tensorMap_N_PlasticData[6]

    plasticVars::PlasticVars = SmallStrainPlastic.initPlasticVars(model)
    StressDim::Int64 = size(C,1)
    âˆ‚Î¾_âˆ‚xFunc::Function = getFunction_âˆ‚Î¾_âˆ‚x(element)
    noOfIpPoints::Int64 = length(shapeFunction)
    noOfNodes::Int64 = size(shapeFunction[1].âˆ‚Ï•_âˆ‚Î¾,1)
    Ïµáµ–_g::Array{Array{Float64,1},1} = Array{Array{Float64,1},1}(undef, noOfIpPoints)
    for ipNo::Int64 âˆˆ 1:noOfIpPoints
        Ïµáµ–_g[ipNo] = zeros(StressDim)
        #âˆ‚x_âˆ‚Î¾::Array{Float64,2} = get_âˆ‚x_âˆ‚Î¾(coordArray, shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾)
        #âˆ‚Î¾_dx::Array{Float64,2} = âˆ‚Î¾_âˆ‚xFunc(âˆ‚x_âˆ‚Î¾)
        #âˆ‚Ï•_âˆ‚x::Array{Float64} = shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾*âˆ‚Î¾_dx
        SmallStrainPlastic.getState!(plasticVars.Ïµáµ–, plasticVars.Î±, stateDict, elementNo, ipNo)
        #println(plasticVars.Ïµáµ–)
        Ïµáµ–_g[ipNo] .= plasticVars.Ïµáµ–
    end
    return Ïµáµ–_g
end

function gaussian_Ïµ(tensorMap_N_PlasticData::T,
    solAtNodes::Array{Float64,1}, problemDim::Int64,
    element::AbstractElement, elementNo::Int64, shapeFunction::Array{ShapeFunction},
    coordArray::Array{Float64,2}; kwargs4function...)::Array{Array{Float64,1},1} where T
    mapDict::Dict{Int64, Int64} = tensorMap_N_PlasticData[1]
    C::Array{Float64,2} = tensorMap_N_PlasticData[2]
    model::PlasticModel = tensorMap_N_PlasticData[3]
    modelParams::ModelParams  = tensorMap_N_PlasticData[4]
    stateDict  = tensorMap_N_PlasticData[5]
    stateDictBuffer  = tensorMap_N_PlasticData[6]

    plasticVars::PlasticVars = SmallStrainPlastic.initPlasticVars(model)
    StressDim::Int64 = size(C,1)
    âˆ‚Î¾_âˆ‚xFunc::Function = getFunction_âˆ‚Î¾_âˆ‚x(element)
    noOfIpPoints::Int64 = length(shapeFunction)
    noOfNodes::Int64 = size(shapeFunction[1].âˆ‚Ï•_âˆ‚Î¾,1)
    Ïµ_g::Array{Array{Float64,1},1} = Array{Array{Float64,1},1}(undef, noOfIpPoints)
    Ïµ::Array{Float64, 1} = zeros(model.ÏµVoigtSize)
    for ipNo::Int64 âˆˆ 1:noOfIpPoints
        Ïµ_g[ipNo] = zeros(StressDim)
        âˆ‚x_âˆ‚Î¾::Array{Float64,2} = get_âˆ‚x_âˆ‚Î¾(coordArray, shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾)
        âˆ‚Î¾_dx::Array{Float64,2} = âˆ‚Î¾_âˆ‚xFunc(âˆ‚x_âˆ‚Î¾)
        âˆ‚Ï•_âˆ‚x::Array{Float64} = shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾*âˆ‚Î¾_dx
        findStrain!(mapDict, Ïµ, âˆ‚Ï•_âˆ‚x,  solAtNodes, problemDim)
        plasticVars.Ïµ = SmallStrainPlastic.getContinuumMandelStrain(Ïµ)
        #println(plasticVars.Ïµáµ–)
        Ïµ_g[ipNo] .= plasticVars.Ïµ
    end
    return Ïµ_g
end

function gaussian_Ïƒ(tensorMap_N_PlasticData::T,
    solAtNodes::Array{Float64,1}, problemDim::Int64,
    element::AbstractElement, elementNo::Int64, shapeFunction::Array{ShapeFunction},
    coordArray::Array{Float64,2}; kwargs4function...)::Array{Array{Float64,1},1} where T
    mapDict::Dict{Int64, Int64} = tensorMap_N_PlasticData[1]
    C::Array{Float64,2} = tensorMap_N_PlasticData[2]
    model::PlasticModel = tensorMap_N_PlasticData[3]
    modelParams::ModelParams  = tensorMap_N_PlasticData[4]
    stateDict  = tensorMap_N_PlasticData[5]
    stateDictBuffer  = tensorMap_N_PlasticData[6]

    plasticVars::PlasticVars = SmallStrainPlastic.initPlasticVars(model)
    plasticVars.C = C
    StressDim::Int64 = size(C,1)
    âˆ‚Î¾_âˆ‚xFunc::Function = getFunction_âˆ‚Î¾_âˆ‚x(element)
    noOfIpPoints::Int64 = length(shapeFunction)
    noOfNodes::Int64 = size(shapeFunction[1].âˆ‚Ï•_âˆ‚Î¾,1)
    Ïƒ_g::Array{Array{Float64,1},1} = Array{Array{Float64,1},1}(undef, noOfIpPoints)
    Ïµ::Array{Float64, 1} = zeros(model.ÏµVoigtSize)
    for ipNo::Int64 âˆˆ 1:noOfIpPoints
        Ïƒ_g[ipNo] = zeros(StressDim)
        âˆ‚x_âˆ‚Î¾::Array{Float64,2} = get_âˆ‚x_âˆ‚Î¾(coordArray, shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾)
        âˆ‚Î¾_dx::Array{Float64,2} = âˆ‚Î¾_âˆ‚xFunc(âˆ‚x_âˆ‚Î¾)
        âˆ‚Ï•_âˆ‚x::Array{Float64} = shapeFunction[ipNo].âˆ‚Ï•_âˆ‚Î¾*âˆ‚Î¾_dx
        findStrain!(mapDict, Ïµ, âˆ‚Ï•_âˆ‚x,  solAtNodes, problemDim)
        plasticVars.Ïµ = SmallStrainPlastic.getContinuumMandelStrain(Ïµ)
        SmallStrainPlastic.getState!(plasticVars.Ïµáµ–, plasticVars.Î±, stateDict, elementNo, ipNo)
        #println(plasticVars.Ïµáµ–)
        Ïƒ_g[ipNo] .= C*(plasticVars.Ïµ - plasticVars.Ïµáµ–)
    end
    return Ïƒ_g
end
#=
j2Model = SmallStrainPlastic.j2Model

function initParams_j2(Ïƒ_y::Float64, params_H::Float64)
    return SmallStrainPlastic.initParams_j2(Ïƒ_y, params_H)
end

function updateStateDict4rmBuffer()
    SmallStrainPlastic.updateStateDict4rmBuffer()
    return nothing
end
=#
