#----------------------------------------------------------------------
#  Direct GFDM method - 2d Navier-Cauchy equations, Mixed BC
#  Uniaxial rod
#  Author: Andrea Pavan
#  Date: 21/12/2022
#  License: GPLv3-or-later
#----------------------------------------------------------------------
using ElasticArrays;
using LinearAlgebra;
using SparseArrays;
using PyPlot;
include("utils.jl");


#problem definition
l1 = 10.0;      #domain x size
l2 = 1.0;       #domain y size
utip = 1e-3;    #forced tip displacement
E = 68e9;       #Young modulus
ν = 0.33;       #Poisson ratio

meshSize = 0.2;     #distance target between internal nodes
surfaceMeshSize = 0.2;      #distance target between boundary nodes
minNeighbors = 8;       #minimum number of neighbors allowed
minSearchRadius = meshSize/2;       #starting search radius


#Euler-Bernoulli solution
println("Euler-Bernoulli beam theory");
println("  utip = ",utip);
println("  σxx = ",E*utip/l1);
μ = 0.5*E/(1+ν);        #Lamè coefficients
λ = E*ν/((1+ν)*(1-2*ν));


#pointcloud generation
time1 = time();
pointcloud = ElasticArray{Float64}(undef,2,0);      #2xN matrix containing the coordinates [X;Y] of each node
boundaryNodes = Vector{Int};        #indices of the boundary nodes
normals = ElasticArray{Float64}(undef,2,0);     #2xN matrix containing the components [nx;ny] of the normal of each boundary node
for i=0:surfaceMeshSize:l1
    append!(pointcloud, [i,0.0]);
    append!(normals, [0,-1]);
end
for i=surfaceMeshSize:surfaceMeshSize:l2
    append!(pointcloud, [l1,i]);
    append!(normals, [1,0]);
end
for i=l1-surfaceMeshSize:-surfaceMeshSize:0
    append!(pointcloud, [i,l2]);
    append!(normals, [0,1]);
end
for i=l2-surfaceMeshSize:-surfaceMeshSize:surfaceMeshSize
    append!(pointcloud, [0,i]);
    append!(normals, [-1,0]);
end
boundaryNodes = collect(range(1,size(pointcloud,2)));
for x=0+meshSize:meshSize:l1-meshSize
    for y=0+meshSize:meshSize:l2-meshSize
        newP = [x,y]+(rand(Float64,2).-0.5).*meshSize/5;
        insertP = true;
        for j in boundaryNodes
            if (newP[1]-pointcloud[1,j])^2+(newP[2]-pointcloud[2,j])^2<(0.75*meshSize)^2
                insertP = false;
            end
        end
        if insertP
            append!(pointcloud, newP);
        end
    end
end
internalNodes = collect(range(1+length(boundaryNodes),size(pointcloud,2)));
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));

#boundary conditions
N = size(pointcloud,2);     #number of nodes
uD = Vector{Float64}(undef,N);
vD = Vector{Float64}(undef,N);
uN = Vector{Float64}(undef,N);
vN = Vector{Float64}(undef,N);
for i=1:N
    uD[i] = NaN;
    vD[i] = NaN;
    uN[i] = NaN;
    vN[i] = NaN;
end
for i in boundaryNodes
    if pointcloud[1,i]==l1
        #right boundary
        uD[i] = utip;
        vN[i] = 0.0;
    else
        if pointcloud[1,i]==0
            #left boundary
            uD[i] = 0.0;
            vD[i] = 0.0;
        else
            #bottom boundary, top boundary
            if pointcloud[2,i]==0 || pointcloud[2,i]==l2
                uN[i] = 0.0;
                vN[i] = 0.0;
            end
        end
    end
end

#=
#pointcloud plot
figure();
plot(pointcloud[1,boundaryNodes],pointcloud[2,boundaryNodes],"r.");
plot(pointcloud[1,internalNodes],pointcloud[2,internalNodes],"k.");
title("Pointcloud plot");
axis("equal");
display(gcf());
=#


#neighbor search
time2 = time();
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
for i=1:N
    searchradius = minSearchRadius;
    while Nneighbors[i]<minNeighbors
        neighbors[i] = Int[];
        #check every other node
        for j=1:N
            if j!=i && all(abs.(pointcloud[:,j]-pointcloud[:,i]).<searchradius)
                push!(neighbors[i],j);
            end
        end
        unique!(neighbors[i]);
        Nneighbors[i] = length(neighbors[i]);
        searchradius += minSearchRadius/2;
    end
end
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
w2 = Vector{Vector{Float64}}(undef,N);      #neighbors weights
for i=1:N
    P[i] = Array{Float64}(undef,2,Nneighbors[i]);
    r2[i] = Vector{Float64}(undef,Nneighbors[i]);
    w2[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]
        w2[i][j] = exp(-1*r2[i][j]/r2max)^2;
    end
end
w2pde = 2.0;        #least squares weight for the pde
w2bc = 2.0;     #least squares weight for the boundary condition


#least square matrix inversion
A = Vector{Matrix}(undef,N);        #least-squares matrices
C = Vector{Matrix}(undef,N);        #derivatives coefficients matrices
for i in internalNodes
    xj = P[i][1,:];
    yj = P[i][2,:];
    V = zeros(Float64,2+2*Nneighbors[i],12);
    for j=1:Nneighbors[i]
        V[j,:] = [1, xj[j], yj[j], xj[j]^2, yj[j]^2, xj[j]*yj[j], 0, 0, 0, 0, 0, 0];
        V[j+Nneighbors[i],:] = [0, 0, 0, 0, 0, 0, 1, xj[j], yj[j], xj[j]^2, yj[j]^2, xj[j]*yj[j]];
    end
    V[1+2*Nneighbors[i],:] = [0, 0, 0, 2*(2+λ/μ), 2, 0, 0, 0, 0, 0, 0, 1+λ/μ];
    V[2+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 1+λ/μ, 0, 0, 0, 2, 2*(2+λ/μ), 0];
    W = Diagonal(vcat(w2[i],w2[i],w2pde,w2pde));
    A[i] = transpose(V)*W*V;
    (Q,R) = qr(A[i]);
    C[i] = inv(R)*transpose(Q)*transpose(V)*W;
    #(L,U) = lu(A[i]);
    #C[i] = inv(U)*inv(L)*transpose(V)*W;
end
for i in boundaryNodes
    xj = P[i][1,:];
    yj = P[i][2,:];
    V = zeros(Float64,4+2*Nneighbors[i],12);
    for j=1:Nneighbors[i]
        V[j,:] = [1, xj[j], yj[j], xj[j]^2, yj[j]^2, xj[j]*yj[j], 0, 0, 0, 0, 0, 0];
        V[j+Nneighbors[i],:] = [0, 0, 0, 0, 0, 0, 1, xj[j], yj[j], xj[j]^2, yj[j]^2, xj[j]*yj[j]];
    end
    V[1+2*Nneighbors[i],:] = [0, 0, 0, 2*(2+λ/μ), 2, 0, 0, 0, 0, 0, 0, 1+λ/μ];
    V[2+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 1+λ/μ, 0, 0, 0, 2, 2*(2+λ/μ), 0];
    if !isnan(uD[i])
        V[3+2*Nneighbors[i],:] = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
        #V = vcat(V,[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]');
    end
    if !isnan(uN[i])
        V[3+2*Nneighbors[i],:] = [0, normals[1,i], normals[2,i], 0, 0, 0, 0, 0, 0, 0, 0, 0];
        #V = vcat(V,[0, normals[1,i], normals[2,i], 0, 0, 0, 0, 0, 0, 0, 0, 0]');
    end
    if !isnan(vD[i])
        V[4+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0];
        #V = vcat(V,[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]');
    end
    if !isnan(vN[i])
        V[4+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 0, 0, normals[1,i], normals[2,i], 0, 0, 0];
        #V = vcat(V,[0, 0, 0, 0, 0, 0, 0, normals[1,i], normals[2,i], 0, 0, 0]');
    end
    W = Diagonal(vcat(w2[i],w2[i],w2pde,w2pde,w2bc,w2bc));
    A[i] = transpose(V)*W*V;
    (Q,R) = qr(A[i]);
    C[i] = inv(R)*transpose(Q)*transpose(V)*W;
    #(L,U) = lu(A[i]);
    #C[i] = inv(U)*inv(L)*transpose(V)*W;
end
println("Inverted least-squares matrices in ", round(time()-time3,digits=2), " s");


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


#linear system solution
time5 = time();
b = zeros(2*N);       #rhs vector
for i in internalNodes
    b[i] = 0;
end
for i in boundaryNodes
    b[i] = 0;
    b[N+i] = 0;
    if !isnan(uD[i])
        b[i] += C[i][1,end-1]*uD[i];
        b[N+i] += C[i][7,end-1]*uD[i];
    end
    if !isnan(uN[i])
        b[i] += C[i][1,end-1]*uN[i];
        b[N+i] += C[i][7,end-1]*uN[i];
    end
    if !isnan(vD[i])
        b[N+i] += C[i][7,end]*vD[i];
        b[i] += C[i][1,end]*vD[i];
    end
    if !isnan(vN[i])
        b[N+i] += C[i][7,end]*vN[i];
        b[i] += C[i][1,end]*vN[i];
    end
end
sol = M\b;
println("Linear system solved in ", round(time()-time5,digits=2), " s");

#displacement plot
u = sol[1:N];
v = sol[N+1:end];
figure();
scatter(pointcloud[1,:],pointcloud[2,:],c=u,cmap="Oranges");
colorbar();
title("Uniaxial rod - x displacement");
axis("equal");
display(gcf());


#von Mises stress
εxx = Vector{Float64}(undef,N);
σxx = Vector{Float64}(undef,N);
for i=1:N
    εxx[i] = 0.0;
    for j=1:Nneighbors[i]
        εxx[i] += C[i][2,j]*u[neighbors[i][j]] + C[i][2,j+Nneighbors[i]]*v[neighbors[i][j]];
    end
end
for i in boundaryNodes
    if !isnan(uD[i])
        εxx[i] += C[i][2,end-1]*uD[i];
    end
    if !isnan(uN[i])
        εxx[i] += C[i][2,end-1]*uN[i];
    end
    if !isnan(vD[i])
        εxx[i] += C[i][2,end]*vD[i];
    end
    if !isnan(vN[i])
        εxx[i] += C[i][2,end]*vN[i];
    end
end
σxx = E*εxx;
println("Stress σxx = ",maximum(σxx));

#stress plot
figure();
scatter(pointcloud[1,:],pointcloud[2,:],c=σxx,cmap="hsv");
colorbar();
title("Uniaxial rod - σxx stress");
axis("equal");
display(gcf());
