include("LinTetra.jl")
include("QTetraInterp.jl")

"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(eF-E(k)).
"""
function QuadTetraÎ˜(Eqtetra,eF,iter)
    FloatType = typeof(float(eF))
    if iter==0
        if maximum(Eqtetra) <= eF
            # qw = FloatType.([-1/120,-1/120,-1/120,-1/120,1/30,1/30,1/30,1/30,1/30,1/30])
            qw = FloatType.(@SArray [0.005208333333333333, 0.005208333333333333, 
            0.005208333333333333,0.005208333333333333, 0.020833333333333332, 
            0.020833333333333332, 0.031249999999999997, 0.031249999999999997, 
            0.020833333333333332, 0.020833333333333332])
        else 
            tmp = @MArray zeros(FloatType,10)
            if minimum(Eqtetra) < eF
                @views for i=1:8
                    res =  LinTetraÎ˜(SVector{4}(Eqtetra[qtetra2tetra[i,:]]),eF)/8
                    for j=1:4
                        @inbounds    tmp[qtetra2tetra[i,j]] += res[j]
                #    @inbounds tmp[qtetra2tetra[i,:]] += LinTetraÎ˜(SVector{4}(Eqtetra[qtetra2tetra[i,:]]) ,eF)/8
                    end
                end
            end
            qw = SArray(tmp)
        end
    else 
        eqtetras = QTetraInterpolation(Eqtetra)
        qweights = @MArray zeros(FloatType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetraÎ˜(eqtetras[i,:],eF,iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end
"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(X(k)).
"""
QuadTetraÎ˜(Xqtetra,iter) = QuadTetraÎ˜(-Xqtetra,zero(Xqtetra[1]),iter)

"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î´(eF-E(k)).
"""
function QuadTetraÎ´(Eqtetra,eF,iter)
    FloatType = typeof(float(eF))
        if iter ==0
            if minimum(Eqtetra) < eF < maximum(Eqtetra)
                tmp = @MArray zeros(FloatType,10)
                @views for i=1:8
                    res =  LinTetraÎ´(SVector{4}(Eqtetra[qtetra2tetra[i,:]]),eF)/8
                    for j=1:4
                        @inbounds    tmp[qtetra2tetra[i,j]] += res[j]
                    # @inbounds  tmp[qtetra2tetra[i,:]] += LinTetraÎ´( SVector{4}(Eqtetra[qtetra2tetra[i,:]]) ,eF)/8
                    end
                end
                qw = SArray(tmp)
            else
                qw = @SArray zeros(FloatType,10)
            end
        else
            eqtetras = QTetraInterpolation(Eqtetra)
            qweights = @MArray zeros(FloatType,8,10)
            for i =1:8
            # for i=1:8    
                qweights[i,:]=QuadTetraÎ´(eqtetras[i,:],eF,iter-1)
            end
            qw = CollectQWeights(SArray(qweights))
        end
    
    # print(size(qw))
    return qw
end
"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î´(X(k)).
"""
QuadTetraÎ´(Xqtetra,iter)=QuadTetraÎ´(-Xqtetra,0.,iter)

"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(eF-E(k))Î´(D(k)).
"""
function QuadTetraÎ˜Î´(Eqtetra,eF,Dqtetra,iter)
    FloatType = typeof(float(eF))
    if iter ==0
        if minimum(Eqtetra) < eF 
            tmp = @MArray zeros(FloatType,10)
            @views for i=1:8
                res = LinTetraÎ˜Î´(SVector{4}(Eqtetra[qtetra2tetra[i,:]]),eF,SVector{4}(Dqtetra[qtetra2tetra[i,:]]))/8
                for j=1:4
                    @inbounds    tmp[qtetra2tetra[i,j]] += res[j]
                end
            end
            qw = SArray(tmp)
        else
            qw = @SArray zeros(FloatType,10)
        end
    else
        eqtetras = QTetraInterpolation(Eqtetra)
        dqtetras = QTetraInterpolation(Dqtetra)
        qweights = @MArray zeros(FloatType,8,10)
        for i =1:8
        # for i=1:8    
            qweights[i,:]=QuadTetraÎ˜Î´(eqtetras[i,:],eF,dqtetras[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end

    return qw
end
"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(X1(k))Î´(X2(k)).
"""
QuadTetraÎ˜Î´(X1qtetra,X2qtetra,iter) = QuadTetraÎ˜Î´(-X1qtetra,zero(X1qtetra[1]),X2qtetra,iter)    


"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(eF-E(k))â‹…1/D(k)
"""
function QuadTetraÎ˜ð�”‡(Eqtetra,eF,Dqtetra,iter)
    FloatType = typeof(float(eF))
    dmax = maximum(Dqtetra)
    dmin = minimum(Dqtetra)
    if dmax*dmin>0 && isapprox(dmax,dmin,rtol = 0.25)
        qw = QuadTetraÎ˜(Eqtetra,eF,iter) ./ Dqtetra
    else 
        if iter ==0
            tmp = @MArray zeros(FloatType,10)
            if minimum(Eqtetra) < eF 
                @views for i=1:8
                    res =  LinTetraÎ˜ð�”‡(SVector{4}(Eqtetra[qtetra2tetra[i,:]]),eF,SVector{4}(Dqtetra[qtetra2tetra[i,:]]))/8
                    for j=1:4
                        @inbounds    tmp[qtetra2tetra[i,j]] += res[j]
                    # LinTetraÎ˜ð�”‡(SVector{4}(Eqtetra[qtetra2tetra[i,:]]),eF,SVector{4}(Dqtetra[qtetra2tetra[i,:]]))/8
                    end
                end
            end
            qw = SArray(tmp)
        else 
            eqtetras = QTetraInterpolation(Eqtetra)
            dqtetras = QTetraInterpolation(Dqtetra)
            qweights = @MArray zeros(FloatType,8,10)
            for i =1:8
                qweights[i,:]=QuadTetraÎ˜ð�”‡(eqtetras[i,:],eF,dqtetras[i,:],iter-1)
            end
            qw = CollectQWeights(SArray(qweights))
        end
    end
    return qw
end
"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(X(k))â‹…1/D(k)
"""
QuadTetraÎ˜ð�”‡(Xqtetra,Dqtetra,iter) = QuadTetraÎ˜ð�”‡(-Xqtetra,zero(Xqtetra[1]),Dqtetra,iter)
