# dependencies
using LFAToolkit
using LinearAlgebra
using Pkg
Pkg.activate("./")
Pkg.instantiate()
using Plots

# setup
p = 3
mesh = Mesh2D(1.0, 1.0)

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# bddc smoother
bddc = LumpedBDDC(diffusion)

# full operator symbols
numbersteps = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2
θ_step = 2π/(numbersteps-1)
θ_range = θ_min:θ_step:θ_max

# compute and plot smoothing factor
# setup
ω = [0.2]
maxeigenvalues = zeros(numbersteps, numbersteps)

# compute
for i in 1:numbersteps, j in 1:numbersteps
    θ = [θ_range[i], θ_range[j]]
    if sqrt(abs(θ[1])^2 + abs(θ[2])^2) > π/128
        A = computesymbols(bddc, ω, θ)
        eigenvalues = [abs(val) for val in eigvals(A)]
        maxeigenvalues[i, j] = max(eigenvalues...)
        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])
    end
end

# plot
println("max eigenvalue: ", maxeigenvalue)
xrange = θ_range/π
heatmap(
    xrange,
    xlabel="θ/π",
    xtickfont=font(12, "Courier"),
    xrange,
    ylabel="θ/π",
    ytickfont=font(12, "Courier"),
    maxeigenvalues,
    title="Lumped BDDC Error Symbol",
    transpose=true,
    aspect_ratio=:equal
)
xlims!(θ_min/π, θ_max/π)
ylims!(θ_min/π, θ_max/π)