#----------------------------------------------------------------------
#  GFDM method test - 2d Navier-Cauchy equation
#  Square plate with circular hole
#  Author: Andrea Pavan
#  Date: 20/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 = 10.0;       #domain y size
a = 1.0;        #hole radius
σ0 = 1.0;       #traction load
E = 1000;       #Young modulus
ν = 0.3;        #Poisson ratio

meshSize = 0.25;
surfaceMeshSize = 0.25;
minNeighbors = 8;
minSearchRadius = meshSize/2;

function exact_σxx(x,y)
    rcoord = sqrt(x^2+y^2);
    θ = atan(y,x);
    return σ0*(1- ((a^2)/(rcoord^2))*(1.5*cos(2*θ)+cos(4*θ))+1.5*((a^4)/(rcoord^4))*cos(4*θ));
end
μ = 0.5*E/(1+ν);      #Lamè coefficients
λ = E*ν/((1+ν)*(1-2*ν));


#pointcloud generation - boundary nodes
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=a+surfaceMeshSize/2:surfaceMeshSize:l1-surfaceMeshSize/2
    append!(pointcloud, [i,0]);
    append!(normals, [0,-1]);
end
for i=0+surfaceMeshSize/2:surfaceMeshSize:l2-surfaceMeshSize/2
    append!(pointcloud, [l1,i]);
    append!(normals, [1,0]);
end
for i=l1-surfaceMeshSize/2:-surfaceMeshSize:0+surfaceMeshSize/2
    append!(pointcloud, [i,l2]);
    append!(normals, [0,1]);
end
for i=l2-surfaceMeshSize/2:-surfaceMeshSize:a+surfaceMeshSize/2
    append!(pointcloud, [0,i]);
    append!(normals, [-1,0]);
end
dθ = surfaceMeshSize/a;
for i=pi/2-dθ/4:-dθ:0+dθ/4
    append!(pointcloud, a.*[cos(i),sin(i)]);
    append!(normals, -[cos(i),sin(i)]);
end
boundaryNodes = collect(range(1,size(pointcloud,2)));

#pointcloud generation - internal nodes (random)
#=NinternalPoints = 0;
surfaceMargin = surfaceMeshSize/2;
while NinternalPoints<1000
    x = rand(-l1/2+surfaceMargin:1e-6:l1/2-surfaceMargin);
    y = rand(0+surfaceMargin:1e-6:l2-surfaceMargin);
    if abs(x)^2+abs(y)^2>(a+surfaceMargin)^2
        append!(pointcloud, [x,y]);
        global NinternalPoints += 1;
    end
end=#

#pointcloud generation - internal nodes (cartesian)
#=surfaceMargin = max(surfaceMeshSize, meshSize);
for x=-l1/2+surfaceMargin:meshSize:l1/2-surfaceMargin
    for y=0+surfaceMargin:meshSize:l2-surfaceMargin
        if abs(x)^2+abs(y)^2>(a+surfaceMargin/2)^2
            #append!(pointcloud, [x,y]);
            append!(pointcloud, [x,y]+(rand(Float64,2).-0.5).*meshSize/5);
        end
    end
end=#

#pointcloud generation - internal nodes (quadtree)
#=surfaceMargin = max(surfaceMeshSize/2, meshSize/2);
(quadtree,quadtreeSize,quadtreeCenter,quadtreePoints,quadtreeNpoints) = buildQuadtree(pointcloud);
quadtreeLeaves = findall(quadtreeNpoints.>=0);
for i in quadtreeLeaves
    #if quadtreeNpoints[i] != 1
        (x,y) = quadtreeCenter[:,i];
        if abs(x)^2+abs(y)^2>=(a+surfaceMargin/2)^2 && x>=-l1/2+surfaceMargin && x<=l1/2-surfaceMargin && y>=0+surfaceMargin && y<=l2-surfaceMargin
            append!(pointcloud, quadtreeCenter[:,i]);
        end
    #end
end=#

#pointcloud generation - internal nodes (layered growth)
#=surfaceMargin = max(surfaceMeshSize/2, meshSize/2);
activeNodes = boundaryNodes;
addedNodes = Vector{Int}(undef,0);
for i=1:length(activeNodes)-1
    B = 0.5.*(pointcloud[:,i]+pointcloud[:,i+1]);
    n = 0.5.*(normals[:,i]+normals[:,i+1]);
    append!(pointcloud, B-n.*surfaceMeshSize./sqrt(2));
    push!(addedNodes,i+length(boundaryNodes));
end
while length(addedNodes)>5
    global activeNodes = copy(addedNodes);
    global addedNodes = Vector{Int}(undef,0);
    for i=1:length(activeNodes)-1
        B = 0.5.*(pointcloud[:,activeNodes[i]]+pointcloud[:,activeNodes[i+1]]);
        n = [pointcloud[2,activeNodes[i]]-pointcloud[2,activeNodes[i+1]], pointcloud[1,activeNodes[i+1]]-pointcloud[1,activeNodes[i]]];
        #n ./= sqrt(n[1]^2+n[2]^2);
        #newP = B+n.*meshSize./sqrt(2);
        newP = B+n./2;
        if newP[1]^2+newP[2]^2>=(a+surfaceMargin/2)^2 && newP[1]>=0+surfaceMargin && newP[1]<=l1+surfaceMargin && newP[2]>=0+surfaceMargin && newP[2]<=l2-surfaceMargin
            insertP = true;
            for j in size(pointcloud,2)
                if (newP[1]-pointcloud[1,j])^2+(newP[2]-pointcloud[2,j])^2<(1.2*meshSize)^2
                    insertP = false;
                    break;
                end
            end
            if insertP
                append!(pointcloud, newP);
                push!(addedNodes,size(pointcloud,2));
            end
        end
    end
    #face [end,1]
    B = 0.5.*(pointcloud[:,activeNodes[end]]+pointcloud[:,activeNodes[1]]);
    n = [pointcloud[2,activeNodes[end]]-pointcloud[2,activeNodes[1]], pointcloud[1,activeNodes[1]]-pointcloud[1,activeNodes[end]]];
    #n ./= sqrt(n[1]^2+n[2]^2);
    #newP = B+n.*meshSize./sqrt(2);
    newP = B+n./2;
    if newP[1]^2+newP[2]^2>=(a+surfaceMargin/2)^2 && newP[1]>=0+surfaceMargin && newP[1]<=l1+surfaceMargin && newP[2]>=0+surfaceMargin && newP[2]<=l2-surfaceMargin
        insertP = true;
        for j in size(pointcloud,2)
            if (newP[1]-pointcloud[1,j])^2+(newP[2]-pointcloud[2,j])^2<(1.2*meshSize)^2
                insertP = false;
                break;
            end
        end
        if insertP
            append!(pointcloud, newP);
            push!(addedNodes,size(pointcloud,2));
        end
    end
    #growing pointcloud plot
    #=println("Added nodes: ",length(addedNodes));
    internalNodes = collect(range(1+length(boundaryNodes),size(pointcloud,2)));
    figure(1);
    plot(pointcloud[1,boundaryNodes],pointcloud[2,boundaryNodes],"r.");
    plot(pointcloud[1,internalNodes],pointcloud[2,internalNodes],"k.");
    title("Pointcloud plot");
    axis("equal");
    display(gcf());
    sleep(2);=#
end=#

#pointcloud generation - internal nodes (flowmesher)
NinternalPoints = 0;
surfaceMargin = surfaceMeshSize/2;
while NinternalPoints<1000
    x = rand(0+surfaceMargin:1e-6:l1-surfaceMargin);
    y = rand(0+surfaceMargin:1e-6:l2-surfaceMargin);
    if x^2+y^2>(a+surfaceMargin)^2
        append!(pointcloud, [x,y]);
        global NinternalPoints += 1;
    end
end
P = pointcloud;

#initial pointcloud plot
figure();
plot(pointcloud[1,boundaryNodes],pointcloud[2,boundaryNodes],"r.");
plot(pointcloud[1,length(boundaryNodes)+1:size(pointcloud,2)],pointcloud[2,length(boundaryNodes)+1:size(pointcloud,2)],"k.");
title("Initial pointcloud plot");
axis("equal");
display(gcf());

#v = rand(-meshSize/5:1e-6:meshSize/5,2,size(P,2));
v = zeros(Float64,2,size(P,2));
F = zeros(Float64,2,size(P,2));
#maxV = maximum(abs.(v));
maxV = 1.0;
dt = 0.1;
h = 3*meshSize;
iter = 1;
while maxV>5e-3
    for i=1+length(boundaryNodes):size(P,2)
        global F[:,i] = [0,0];
        #if v[:,i]'v[:,i]<1e-2
            for j=1:size(P,2)
                rd2 = transpose(P[:,j]-P[:,i])*(P[:,j]-P[:,i]);
                if i!=j && rd2<h^2
                    global F[:,i] += exp(-6.5*rd2/(h^2))*(P[:,i]-P[:,j])/sqrt(rd2);
                end
            end
        #end
        ang = atan(P[2,i],P[1,i]);
        #F[:,i] .+= 1e-1/(transpose(r1.*[cos(ang),sin(ang)]-P[:,i])*(r1.*[cos(ang),sin(ang)]-P[:,i]));
        #F[:,i] .-= 1e-1./(transpose(r2.*[cos(ang),sin(ang)]-P[:,i])*(r2.*[cos(ang),sin(ang)]-P[:,i]));
    end
    global v += dt*F;
    for i=1+length(boundaryNodes):size(P,2)
        modV = sqrt(v[1,i]^2+v[2,i]^2);
        if modV>0.5*meshSize/dt
            v[:,i] /= modV;
            v[:,i] *= 0.5*meshSize/dt;
        end
    end
    global v .*= 0.5;
    global maxV = sum(abs.(v[:]))/length(v[:]);
    println("Flowmesher: iter = ", iter, ", meanV = ",round(maxV,digits=4));
    global P += dt.*v;
    for i=1+length(boundaryNodes):size(P,2)
        if P[1,i]^2+P[2,i]^2<(a+surfaceMargin)^2
            ang = atan(P[2,i],P[1,i]);
            rad = sqrt(P[:,i]'P[:,i]);
            P[:,i] += 2*((a+surfaceMargin).*[cos(ang),sin(ang)]-P[:,i]);
            v[:,i] = [0,0];
        end
        if P[1,i]>l1-surfaceMargin
            P[1,i] -= 2*(P[1,i]-l1+surfaceMargin);
            v[:,i] = [0,0];
        end
        if P[1,i]<0+surfaceMargin
            P[1,i] -= 2*(P[1,i]-0-surfaceMargin);
            v[:,i] = [0,0];
        end
        if P[2,i]>l2-surfaceMargin
            P[2,i] -= 2*(P[2,i]-l2+surfaceMargin);
            v[:,i] = [0,0];
        end
        if P[2,i]<0+surfaceMargin
            P[2,i] -= 2*(P[2,i]-0-surfaceMargin);
            v[:,i] = [0,0];
        end
    end
    global iter += 1;
    #=if mod(iter,10)==0
        figure();
        #plot(pointcloud[1,boundaryNodes],pointcloud[2,boundaryNodes],"r.");
        plot(P[1,:],P[2,:],"k.");
        title("Pointcloud plot");
        axis("equal");
        display(gcf());
        sleep(0.1);
    end=#
end
pointcloud = P;

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));

#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());


#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;
    uN[i] = 0.0;
    vD[i] = NaN;
    vN[i] = 0.0;
end
for i in boundaryNodes
    if pointcloud[2,i]==0
        #bottom
        vD[i] = 0.0;
    end
    if pointcloud[1,i]==l1
        #right
        uN[i] = σ0;
    end
    if pointcloud[1,i]==0
        #left
        uD[i] = 0.0;
    end
end


#neighbor search
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
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),")");

#connectivity plot
#=figure();
plot2Idx = rand(1:N,5);
plot(pointcloud[1,:],pointcloud[2,:],marker=".",linestyle="None",color="lightgray");
for i in plot2Idx
    connColor = rand(3);
    plot(pointcloud[1,neighbors[i]],pointcloud[2,neighbors[i]],marker=".",linestyle="None",color=connColor);
    for j in neighbors[i]
        plot([pointcloud[1,i],pointcloud[1,j]],[pointcloud[2,i],pointcloud[2,j]],"-",color=connColor);
    end
end
plot(pointcloud[1,plot2Idx],pointcloud[2,plot2Idx],"k.");
title("Connectivity plot");
axis("equal");
display(gcf());=#

#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(-6*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;
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];
    else
        #V[3+2*Nneighbors[i],:] = [0, normals[1,i], normals[2,i], 0, 0, 0, 0, 0, 0, 0, 0, 0];
        #V[3+2*Nneighbors[i],:] = [0, normals[1,i]*(λ/μ)*(1-ν)/ν, 2*normals[2,i], 0, 0, 0, 0, 0, normals[1,i]*(λ/μ), 0, 0, 0];
        V[3+2*Nneighbors[i],:] = [0, normals[1,i]*(2+λ/μ), normals[2,i], 0, 0, 0, 0, normals[2,i], normals[1,i]*(λ/μ), 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];
    else
        #V[4+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 0, 0, normals[1,i], normals[2,i], 0, 0, 0];
        #V[4+2*Nneighbors[i],:] = [0, normals[2,i]*(λ/μ), 0, 0, 0, 0, 0, 2*normals[1,i], normals[2,i]*(λ/μ)*(1-ν)/ν, 0, 0, 0];
        V[4+2*Nneighbors[i],:] = [0, normals[2,i]*(λ/μ), normals[1,i], 0, 0, 0, 0, normals[1,i], normals[2,i]*(2+λ/μ), 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;
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];
    else
        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];
    else
        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("Hole plate - x displacement");
axis("equal");
display(gcf());=#


#σxx stress
dudx = Vector{Float64}(undef,N);
dvdy = Vector{Float64}(undef,N);
σxx = Vector{Float64}(undef,N);
for i=1:N
    dudx[i] = 0.0;
    for j=1:Nneighbors[i]
        dudx[i] += C[i][2,j]*u[neighbors[i][j]] + C[i][2,j+Nneighbors[i]]*v[neighbors[i][j]];
        dvdy[i] += C[i][9,j]*u[neighbors[i][j]] + C[i][9,j+Nneighbors[i]]*v[neighbors[i][j]];
    end
end
for i in boundaryNodes
    if !isnan(uD[i])
        dudx[i] += C[i][2,end-1]*uD[i];
        dvdy[i] += C[i][9,end-1]*uD[i];
    else
        dudx[i] += C[i][2,end-1]*uN[i]/μ;
        dvdy[i] += C[i][9,end-1]*uN[i]/μ;
    end
    if !isnan(vD[i])
        dudx[i] += C[i][2,end]*vD[i];
        dvdy[i] += C[i][9,end]*vD[i];
    else
        dudx[i] += C[i][2,end]*vN[i]/μ;
        dvdy[i] += C[i][9,end]*vN[i]/μ;
    end
end
#σxx = dudx*(1-ν)*E/((1+ν)*(1-2ν)) + dvdy*ν*E/((1+ν)*(1-2ν));
#σxx = E*dudx;
#σxx = dudx*λ*(1-ν)/ν;
#σxx = dudx*λ*(1-ν)/ν + dvdy*λ;
σxx = dudx*(2μ+λ) + dvdy*λ;

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


#validation plot - σxx(0,y) stress
idxPlot = findall(pointcloud[1,:].==0);
y_exactsol = collect(a:0.01:l2);
exactsol = 0*y_exactsol;
for i=1:length(exactsol)
    exactsol[i] = exact_σxx(0,y_exactsol[i]);
end
σxx_rmse = sqrt(sum((σxx-exact_σxx.(pointcloud[1,:],pointcloud[2,:])).^2)/N);
figure();
plot(pointcloud[2,idxPlot],σxx[idxPlot],"r.",label="GFDM");
plot(y_exactsol,exactsol,"k-",linewidth=1.0,label="Analytical");
title("σxx stress @x=0");
legend(loc="upper right");
xlabel("y coordinate");
ylabel("σxx stress");
axis("equal");
display(gcf());

#validation plot - σxx(x,0) stress
idxPlot = findall(pointcloud[2,:].==0);
y_exactsol = collect(a:0.01:l1);
exactsol = 0*y_exactsol;
for i=1:length(exactsol)
    exactsol[i] = exact_σxx(y_exactsol[i],0);
end
figure();
plot(pointcloud[1,idxPlot],σxx[idxPlot],"r.",label="GFDM");
plot(y_exactsol,exactsol,"k-",linewidth=1.0,label="Analytical");
title("σxx stress @y=0");
legend(loc="upper right");
xlabel("x coordinate");
ylabel("σxx stress");
axis("equal");
display(gcf());

#exact solution plot
#=figure();
scatter(pointcloud[1,:],pointcloud[2,:],c=exact_σxx.(pointcloud[1,:],pointcloud[2,:]),cmap="jet");
colorbar();
title("Exact solution");
axis("equal");
display(gcf());=#
