#====================================================================
  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/.
 =====================================================================#




#=This algorithm is as per the algorithm mentioned in The Finite Element
Method: Its Basis and Fundamentals, Sixth edition, O.C. Zienkiewicz
Ch-17 The time dimension â€“ discrete approximation in time=#

"""As per the mentioned reference, the A matrix is calcualated as,

A = Î¸[p-2]*Î”t^(p-2)/factorial(p-2)*M + Î¸[p-1]*Î”t^(p-1)/factorial(p-1)*K + Î¸[p]*Î”t^(p)/factorial(p)*C

There is a difference between the reference formula in Zienkiewicz's book and the one used since a Î¸[0] = 1 is not
possible in julia and Î¸[1] = 1 is taken rest of elements are the same as Î¸_Array.
"""
function build_A(M_C_K_Array::Array{SparseMatrixCSC{Float64,Int64},1}, Î”t::Float64, thetaModified::Array{Float64, 1})
    @assert length(M_C_K_Array) <= length(thetaModified) "Length of Î¸_array must be less 1 than, equal to, or greater than M_C_K_Array"
    A::SparseMatrixCSC = spzeros(size(M_C_K_Array[1])...)

    p::Int64 = length(thetaModified) - length(M_C_K_Array) + 1
    for n âˆˆ 1:length(M_C_K_Array)
        A += (thetaModified[p]*Î”t^(p-1)/float(factorial(p-1)))*M_C_K_Array[n]
        p += 1
    end
    return A
end

"""This function evaluates the values of Ì‡fÌ‡, fÌˆ, ..., etc. For this it uses taylor series and the values of fâ‚™, fâ‚™â‚Šâ‚�, fâ‚™â‚Šâ‚‚, etc."""
function taylor_f_Dot_eval(f_Array::Array{Array{Float64,1},1}, Î”t::Float64)
    rowlength::Int64 = length(f_Array)-1
    f_n_Matrix::Array{Float64, 2} = Array{Float64, 2}(undef, rowlength, length(f_Array[1]))
    for i âˆˆ 1:rowlength
        f_n_Matrix[i,:] = (f_Array[i+1]-f_Array[1])'
    end
    Î”t_Matrix::Array{Float64, 2} = Array{Float64, 2}(undef, rowlength, rowlength)
    Î”t_Matrix[:,1] = Î”t .*collect(1:rowlength)
    for i âˆˆ 2:rowlength
        Î”t_Matrix[:,i] = (Î”t_Matrix[:,1].^i)/float(factorial(i))
    end
    f_DotMatrix::Array{Float64, 2} = Î”t_Matrix\f_n_Matrix
    #println("f_DotMatrix = ",f_DotMatrix)
    return f_DotMatrix
end

"""This function evaluates the mean values of Ì‡the solution array, u_mean, uÌ‡_mean, uÌˆ_mean, etc """
function buildSolMeanArray(SolutionArray::Array{Array{Float64, 1},1}, Î”t::Float64,
    Î¸_Array::Array{Float64,1})

    SolMeanArray::Array{Array{Float64, 1},1} = Array{Array{Float64, 1},1}(undef, length(SolutionArray))
    for dot âˆˆ 1:length(SolutionArray)
        p::Int64 = 1
        SolMeanArray[dot] = zeros(length(SolutionArray[1]))
        SolMeanArray[dot] = SolutionArray[dot]
        for i âˆˆ dot+1:length(SolutionArray)
            #println(p)
            SolMeanArray[dot] += Î¸_Array[p]*Î”t^(p)/factorial(p)*SolutionArray[i]
            p +=1
        end
        #p = 0
    end
    #println("SolMeanArray = ", SolMeanArray)
    return SolMeanArray
end

"""Builds f_mean so the solution can be evaluated as Î± = Aâ�»Â¹f_mean"""
function build_f_mean(f_Array::Array{Array{Float64,1},1}, Î”t::Float64, Î¸_Array::Array{Float64, 1})
    @assert length(f_Array) == length(Î¸_Array)+1 "Length of f_Array must be 1 greater than length of Î¸_Array"
    f_DotMatrix::Array{Float64, 2} = taylor_f_Dot_eval(f_Array, Î”t)
    f_mean::Array{Float64, 1} = deepcopy(f_Array[1])
    factorial_i::Array{Float64, 1} = Array{Float64, 1}(undef, size(f_DotMatrix,1))
    for i âˆˆ 1:size(f_DotMatrix,1)
        factorial_i[i] = float(factorial(i))
    end
    for j âˆˆ 1:size(f_DotMatrix,2)
        for i âˆˆ 1:size(f_DotMatrix,1)
            f_mean[j] += Î¸_Array[i]*Î”t^(i)/factorial_i[i]*f_DotMatrix[i,j]
        end
    end
    #println("f_mean = ", f_mean)
    return f_mean
end

"""Return A matrix and f_mean vector so that the form ð��€Î± + ð��ŒuÌˆ_mean + ð��‚uÌ‡_mean + ð��Šu_mean + ð��Ÿ_mean = 0 can be formed.

    A::SparseMatrixCSC, SolMeanArray::Array{Array{Float64, 1},1}, f_mean::Array{Float64, 1} = get_SSpj_A_u_f(M_C_K_Array, SolutionArray, f_Array, Î”t, Î¸_Array)
"""
function get_SSpj_A_meanU_f(M_C_K_Array::Array{SparseMatrixCSC{Float64,Int64},1}, SolutionArray::Array{Array{Float64, 1},1},
    f_Array::Array{Array{Float64,1},1}, Î”t::Float64, Î¸_Array::Array{Float64,1})

    thetaModified::Array{Float64, 1} = Array{Float64, 1}(undef, length(Î¸_Array)+1)
    thetaModified[1] = 1.0
    thetaModified[2:end] = Î¸_Array
    A::SparseMatrixCSC = build_A(M_C_K_Array, Î”t, thetaModified)
    f_mean::Array{Float64, 1} = build_f_mean(f_Array, Î”t, Î¸_Array)
    SolMeanArray::Array{Array{Float64, 1},1} = buildSolMeanArray(SolutionArray, Î”t, Î¸_Array)
    return A, SolMeanArray, f_mean
end

"""Returns Sparse Matrix A and vector b so that Î± = Aâ�»Â¹b can be evaluated.

    A, b = SSpj_getFinal_A_b(M_C_K_Array, SolutionArray, f_Array, Î”t, Î¸_Array)
"""
function SSpj_getFinal_A_b(M_C_K_Array::Array{SparseMatrixCSC{Float64,Int64},1}, SolutionArray::Array{Array{Float64, 1},1},
    f_Array::Array{Array{Float64,1},1}, Î”t::Float64, Î¸_Array::Array{Float64,1})

    A::SparseMatrixCSC, SolMeanArray::Array{Array{Float64, 1},1}, f_mean::Array{Float64, 1} = get_SSpj_A_meanU_f(
    M_C_K_Array, SolutionArray, f_Array, Î”t, Î¸_Array)
    lnM_C_K_Array::Int64 = length(M_C_K_Array)
    lnSolMeanArray::Int64 = length(SolutionArray)
    for i âˆˆ 1:lnSolMeanArray
        f_mean += M_C_K_Array[lnM_C_K_Array]*SolMeanArray[i]
        #println("f_mean at i: ", i, " =  ", f_mean[1])
        lnM_C_K_Array -=1
    end
    return -A, f_mean
end

"""Updates solutions: u, uÌ‡, uÌˆ, etc; after obtaining the solution of Î±

    updateSolution!(SolutionArray::Array{Array{Float64, 1},1}, Î”t::Float64, Î±::Array{Float64, 1})
"""
function updateSolution!(SolutionArray::Array{Array{Float64, 1},1}, Î”t::Float64, Î±::Array{Float64, 1})

    for dot âˆˆ 1:length(SolutionArray)
        p::Int64 = 1
        for i âˆˆ dot+1:length(SolutionArray)#(dot-1)
            SolutionArray[dot] += Î”t^p/factorial(p)*SolutionArray[i]
            p +=1
        end
        SolutionArray[dot] += Î”t^p/factorial(p)*Î±
    end
    return nothing
end


"""Updates array of f vectors after obtaining the solution of f at a new time step.

    update_f!(f_Array::Array{Array{Float64,1},1}, f_new::Array{Float64,1})
"""
function update_f!(f_Array::Array{Array{Float64,1},1}, f_new::Array{Float64,1})
    popfirst!(f_Array)
    push!(f_new)
    return nothing
end
