#----------------------------------------------------------------------
#  GFDM method test - 3d oblong plastic cooling
#  Author: Andrea Pavan
#  Date: 06/01/2023
#  License: GPLv3-or-later
#----------------------------------------------------------------------
using ElasticArrays;
using LinearAlgebra;
using SparseArrays;
using PyPlot;
include("utils.jl");


#problem definition
l1 = 0.55e-3;       #x domain size (extrusion width)
l2 = 0.2e-3;        #y domain size (layer height)
l3 = 1e-3;      #z domain size (extrusion length)
kconst = 0.24;      #thermal conductivity
ρ = 1370;       #material density
cp = 1050;      #specific heat capacity
Tamb = 20;      #ambient temperature
Tbed = 50;      #bed temperature
hconst = 20;        #convective film coefficient
T0 = 320;       #initial temperature
t0 = 0;     #starting time
Δt = 0.05;      #timestep
tf = 0.5;       #ending time

meshSize = 0.05e-3;
#meshSize = 0.02e-3;
surfaceMeshSize = meshSize;
minNeighbors = 30;
minSearchRadius = meshSize;


#read pointcloud from a SU2 file
time1 = time();
pointcloud = ElasticArray{Float64}(undef,3,0);      #3xN matrix containing the coordinates [X;Y;Z] of each node
boundaryNodes = Vector{Int}(undef,0);       #indices of the boundary nodes
internalNodes = Vector{Int}(undef,0);       #indices of the internal nodes
normals = ElasticArray{Float64}(undef,3,0);     #3xN matrix containing the components [nx;ny;nz] of the normal of each boundary node

pointcloud = parseSU2mesh("17_3d_heat_filament_1613.su2");
#pointcloud = parseSU2mesh("17_3d_heat_filament_14969.su2");
#cornerPoint = findall((pointcloud[2,:].<=1e-6).*(pointcloud[1,:].<=1e-6));
#pointcloud = pointcloud[:, setdiff(1:end,cornerPoint)];
N = size(pointcloud,2);
idxC = 0;       #index of the reference point
for i=1:N
    if pointcloud[2,i]>=l2/2-1e-6
        #top surface
        push!(boundaryNodes, i);
        append!(normals, [0,1,0]);
        if 0.5*(l1-l2)-1e-6<=pointcloud[1,i]<=0.5*(l1-l2)+1e-6 && pointcloud[3,i]<=0+1e-6
            global idxC = i;
            println("Reference point C index: ",idxC);
        end
    elseif pointcloud[2,i]<=-l2/2+1e-6
        #bottom surface
        push!(boundaryNodes, i);
        append!(normals, [0,-1,0]);
    elseif pointcloud[1,i]>=(l1-l2)/2-1e-6 && (pointcloud[1,i]-(l1-l2)/2)^2+pointcloud[2,i]^2>=(l2/2)^2
        #right bulb
        push!(boundaryNodes, i);
        append!(normals, [pointcloud[1,i]-(l1-l2)/2,pointcloud[2,i],0]./(l2/2));
    elseif pointcloud[1,i]<=-(l1-l2)/2+1e-6 && (pointcloud[1,i]+(l1-l2)/2)^2+pointcloud[2,i]^2>=(l2/2)^2
        #left bulb
        push!(boundaryNodes, i);
        append!(normals, [pointcloud[1,i]+(l1-l2)/2,pointcloud[2,i],0]./(l2/2));
    elseif pointcloud[3,i]<=1e-6
        #rear section
        push!(boundaryNodes, i);
        append!(normals, [0,0,-1]);
    elseif pointcloud[3,i]>=l3-1e-6
        #front section
        push!(boundaryNodes, i);
        append!(normals, [0,0,1]);
    else
        push!(internalNodes, i);
        append!(normals, [0,0,0]);
    end
end

println("Generated pointcloud in ", round(time()-time1,digits=2), " s");
println("Pointcloud properties:");
println("  Boundary nodes: ",length(boundaryNodes));
println("  Internal nodes: ",length(internalNodes));
println("  Memory: ",memoryUsage(pointcloud,boundaryNodes));

#normals plot
#=import Plots;
plt1Idx = boundaryNodes;
figure();
plt1 = Plots.quiver(pointcloud[1,plt1Idx],pointcloud[2,plt1Idx],pointcloud[3,plt1Idx],
        quiver = (normals[1,plt1Idx]/20000,normals[2,plt1Idx]/20000,normals[3,plt1Idx]/20000),
        title = "Normals plot",
        xlim = (-l1,l1),
        ylim = (-l1,l1),
        zlim = (0,l3));
display(plt1);=#


#boundary conditions
N = size(pointcloud,2);     #number of nodes
g1 = zeros(Float64,N);
g2 = zeros(Float64,N);
g3 = zeros(Float64,N);
for i in boundaryNodes
    if pointcloud[2,i]<=-l2/2+1e-6
        #bottom surface
        g1[i] = 1.0;
        g2[i] = 0.0;
        g3[i] = Tbed;
    else
        #everywhere else
        g1[i] = hconst;
        g2[i] = kconst;
        g3[i] = hconst*Tamb;
    end
end

#boundary conditions plot
#=figure();
plt = scatter3D(pointcloud[1,:],pointcloud[2,:],pointcloud[3,:],c=g1);
title("Boundary conditions plot");
axis("equal");
colorbar(plt);
display(gcf());=#


#neighbor search (cartesian cells)
time2 = time();
N = size(pointcloud,2);     #number of nodes
neighbors = Vector{Vector{Int}}(undef,N);       #vector containing N vectors of the indices of each node neighbors
Nneighbors = zeros(Int,N);      #number of neighbors of each node
(neighbors,Nneighbors,cell) = cartesianNeighborSearch(pointcloud,meshSize,minNeighbors);

println("Found neighbors in ", round(time()-time2,digits=2), " s");
println("Connectivity properties:");
println("  Max neighbors: ",maximum(Nneighbors)," (at index ",findfirst(isequal(maximum(Nneighbors)),Nneighbors),")");
println("  Avg neighbors: ",round(sum(Nneighbors)/length(Nneighbors),digits=2));
println("  Min neighbors: ",minimum(Nneighbors)," (at index ",findfirst(isequal(minimum(Nneighbors)),Nneighbors),")");


#neighbors distances and weights
time3 = time();
P = Vector{Array{Float64}}(undef,N);        #relative positions of the neighbors
r2 = Vector{Vector{Float64}}(undef,N);      #relative distances of the neighbors
w = Vector{Vector{Float64}}(undef,N);       #neighbors weights
for i=1:N
    P[i] = Array{Float64}(undef,3,Nneighbors[i]);
    r2[i] = Vector{Float64}(undef,Nneighbors[i]);
    w[i] = Vector{Float64}(undef,Nneighbors[i]);
    for j=1:Nneighbors[i]
        P[i][:,j] = pointcloud[:,neighbors[i][j]]-pointcloud[:,i];
        r2[i][j] = P[i][:,j]'P[i][:,j];
    end
    r2max = maximum(r2[i]);
    for j=1:Nneighbors[i]
        w[i][j] = exp(-6*r2[i][j]/r2max);
        #w[i][j] = 1.0;
    end
end
wpde = 2.0;       #least squares weight for the pde
wbc = 2.0;        #least squares weight for the boundary condition


#least square matrix inversion
C = Vector{Matrix}(undef,N);        #stencil coefficients matrices
#condV = zeros(N);
#condC = zeros(N);
for i in internalNodes
    xj = P[i][1,:];
    yj = P[i][2,:];
    zj = P[i][3,:];
    V = zeros(Float64,1+Nneighbors[i],10);
    for j=1:Nneighbors[i]
        V[j,:] = [1, xj[j], yj[j], zj[j], xj[j]^2, yj[j]^2, zj[j]^2, xj[j]*yj[j], xj[j]*zj[j], yj[j]*zj[j]];
    end
    V[1+Nneighbors[i],:] = [0, 0, 0, 0, 2*kconst/(ρ*cp), 2*kconst/(ρ*cp), 2*kconst/(ρ*cp), 0, 0, 0];
    W = Diagonal(vcat(w[i],wpde));
    #condV[i] = cond(V);
    #(Q,R) = qr(W*V);
    #C[i] = inv(R)*transpose(Matrix(Q))*W;
    VF = svd(W*V);
    C[i] = transpose(VF.Vt)*inv(Diagonal(VF.S))*transpose(VF.U)*W;
    #C[i] = transpose(VF.Vt)*inv(Diagonal(VF.S.+1e-12))*transpose(VF.U)*W;
    #condC[i] = cond(C[i]);
end
for i in boundaryNodes
    #println("Boundary node: ",i);
    xj = P[i][1,:];
    yj = P[i][2,:];
    zj = P[i][3,:];
    V = zeros(Float64,2+Nneighbors[i],10);
    for j=1:Nneighbors[i]
        V[j,:] = [1, xj[j], yj[j], zj[j], xj[j]^2, yj[j]^2, zj[j]^2, xj[j]*yj[j], xj[j]*zj[j], yj[j]*zj[j]];
    end
    V[1+Nneighbors[i],:] = [0, 0, 0, 0, 2*kconst/(ρ*cp), 2*kconst/(ρ*cp), 2*kconst/(ρ*cp), 0, 0, 0];
    V[2+Nneighbors[i],:] = [g1[i], g2[i]*normals[1,i], g2[i]*normals[2,i], g2[i]*normals[3,i], 0, 0, 0, 0, 0, 0];
    W = Diagonal(vcat(w[i],wpde,wbc));
    #condV[i] = cond(V);
    #(Q,R) = qr(W*V);
    #C[i] = inv(R)*transpose(Matrix(Q))*W;
    VF = svd(W*V);
    C[i] = transpose(VF.Vt)*inv(Diagonal(VF.S))*transpose(VF.U)*W;
    #C[i] = transpose(VF.Vt)*inv(Diagonal(VF.S.+1e-12))*transpose(VF.U)*W;
    #condC[i] = cond(C[i]);
end
println("Inverted least-squares matrices in ", round(time()-time3,digits=2), " s");
#=println("Matrix properties:");
println("  Max cond(V): ",maximum(condV));
println("  Avg cond(V): ",sum(condV)/N);
println("  Min cond(V): ",minimum(condV));
println("  Max cond(C): ",maximum(condC));
println("  Avg cond(C): ",sum(condC)/N);
println("  Min cond(C): ",minimum(condC));=#

#condition number distribution
#=figure();
plt.hist(log10.(condC),10);
title("Condition number distribution");
xlabel("log10(cond(C))");
ylabel("Absolute frequency");
display(gcf());=#


#matrix assembly
time4 = time();
rows = Int[];
cols = Int[];
vals = Float64[];
for i=1:N
    push!(rows, i);
    push!(cols, i);
    push!(vals, 1-C[i][1,1+Nneighbors[i]]/Δt);
    for j=1:Nneighbors[i]
        push!(rows, i);
        push!(cols, neighbors[i][j]);
        push!(vals, -C[i][1,j]);
    end
end
M = sparse(rows,cols,vals,N,N);
println("Completed matrix assembly in ", round(time()-time4,digits=2), " s");


#time propagation
time5 = time();
uprev = T0*ones(N);     #solution at the previous step
uC = [];
for t=t0:Δt:tf
    b = zeros(N);       #rhs vector
    for i in internalNodes
        b[i] = -C[i][1,1+Nneighbors[i]]*uprev[i]/Δt;
    end
    for i in boundaryNodes
        b[i] = -C[i][1,1+Nneighbors[i]]*uprev[i]/Δt + C[i][1,2+Nneighbors[i]]*g3[i];
    end
    u = M\b;
    push!(uC, uprev[idxC]);
    println("t = ",t,"; uC = ",uC[end]);
    global uprev = u;
end
println("Simulation completed in ", round(time()-time5,digits=2), " s");

#solution plot
figure();
#plt = scatter3D(pointcloud[1,internalNodes],pointcloud[2,internalNodes],pointcloud[3,internalNodes],c=uprev[internalNodes],cmap="inferno");
plt = scatter3D(pointcloud[1,:],pointcloud[2,:],pointcloud[3,:],c=uprev,cmap="jet");
title("Numerical solution");
axis("equal");
colorbar(plt);
display(gcf());

#validation plot
figure();
plot(collect(t0:Δt:tf),uC,"r-",linewidth=1.0,label="GFDM");
plot([0.01, 0.02, 0.035, 0.0575, 0.09125, 0.14125, 0.19125, 0.24125, 0.29125, 0.34125, 0.39125, 0.44125, 0.49125, 0.5],[314.9072, 306.2313, 286.3006, 252.5652, 207.8958, 159.6279, 125.6891, 102.0784, 85.70526, 74.36204, 66.50572, 61.06476, 57.29651, 56.68461],"k.",label="FEM");
title("Temperature evolution over time");
legend(loc="upper right");
xlabel("Time t");
ylabel("Temperature TC");
display(gcf());
