#----------------------------------------------------------------------
#  GFDM simulation - 3d printing of a narrow wall
#  Author: Andrea Pavan
#  Date: 06/03/2023
#  License: GPLv3-or-later
#----------------------------------------------------------------------
using ElasticArrays;
using LinearAlgebra;
using SparseArrays;
using PyPlot;
include("utils.jl");
include("parseGcode.jl");


#import gcode toolpath
M = parseGcodeFromFile("23_printing_2d_wall_10_008_0055.gcode");
#M = parseGcodeFromFile("23_printing_2d_wall_10_05.gcode");
#M = parseGcodeFromFile("23_printing_2d_wall_70_20.gcode");
#M[:,1] .-= 0.1348;
(lineLength,t,lineFlow,layers,layerHeight) = calculateGcodeFlow(M);
println("Print time: ",round(t[end],digits=0)," s");
println("Filament used: ",round(maximum(M[:,5]),digits=4)," m");

#3d printing process parameters
filamentDiameter = 2.85e-3;
Nsteps = length(lineLength);
lineWidth = 0.4e-3*ones(Nsteps);


#toolpath plot
#=figure();
plot3D(M[:,1],M[:,2],M[:,3]);
title("GCODE tool path");
axis("equal");
display(gcf());=#

#extrusion flow history
#=figure();
plot(t,lineFlow);
title("Extrusion flow history Q(t)");
xlabel("Time t (s)");
xlabel("Extrusion flow (mm/s)");
display(gcf());=#


#simulation
time0 = time();
resolution = 1.0e-3;
meshSize = 0.14e-3/2;
#resolution = 1.0e-3;
#meshSize = 0.15e-3/2;
minNeighbors = 8;
minSearchRadius = 1.1*meshSize;
savePlots = false;
if savePlots
    mkpath("vid1");
    mkpath("vid2");
end

k = 0.24;           #thermal conductivity
ρ = 1370;           #material density
cp = 1050;          #specific heat capacity
Tamb = 20;          #ambient temperature
#Tbed = 50;          #bed temperature
Tbed = 70;
h = 20;             #convective film coefficient
T0 = 230;           #initial temperature
E = 3.1e9;          #Young modulus
ν = 0.43;           #Poisson's ratio
α = 60e-5;          #thermal expansion coefficient
Tref = 230;          #thermal expansion reference temperature

#initialize variables
pointcloud = ElasticArray{Float64}(undef,2,0);      #2xN matrix with [X;Y;Z] as rows
normals = ElasticArray{Float64}(undef,2,0);         #2xN matrix containing the components [nx;ny] of the normal of each boundary node
boundaryNodes = Int[];                              #indices of the boundary nodes
internalNodes = Int[];                              #indices of the internal nodes
depositionTime = Float64[];                         #deposition time for each node
neighbors = Vector{Int}[];                          #vector containing N vectors of the indices of each node neighbors
Nneighbors = Int[];                                 #number of neighbors of each node
N = 0;                                              #number of nodes
Nprev = 0;                                          #number of nodes at the previous timestep
r2max = Float64[];                                  #neighbors star radius
w = Vector{Vector{Float64}}(undef,0);               #neighbors weights
wpde = 2.0;                                         #least squares weight for the governing equation
wbc = 2.0;                                          #least squares weight for the boundary conditions
g1 = Float64[];                                     #boundary conditions
g2 = Float64[];
g3 = Float64[];
g1u = Float64[];
g2u = Float64[];
g3u = Float64[];
g1v = Float64[];
g2v = Float64[];
g3v = Float64[];
dTdx = zeros(Float64,N);
dTdy = zeros(Float64,N);
qx = zeros(Float64,N);
qy = zeros(Float64,N);
C = Vector{Matrix}(undef,0);                        #stencil coefficients matrices
Cm = Vector{Matrix}(undef,0);
Δt = t[2]-t[1];                                     #timestep
T = T0*ones(N);                                     #solution at the previous step
u = zeros(Float64,N);
v = zeros(Float64,N);
dudx = zeros(Float64,N);
dudy = zeros(Float64,N);
dvdx = zeros(Float64,N);
dvdy = zeros(Float64,N);
σxx = zeros(Float64,N);
σyy = zeros(Float64,N);
τxy = zeros(Float64,N);
σv = zeros(Float64,N);
iter = 0;                                           #step number

#time propagation
for τ=2:Nsteps
#for τ=2:30
    if layers[τ]>layers[τ-1]
        println("Layer ",layers[τ], "/",maximum(layers)," (Z = ",round(M[τ,3]*1000,digits=1)," mm)");
    end
    if lineFlow[τ]>1e-9 && lineLength[τ]>1e-9
        svec = (M[τ,1:3]-M[τ-1,1:3])./lineLength[τ];                #direction vector
        Nsubsectors = max(1,round(lineLength[τ]/resolution));       #divide the τ-th step into subsectors depending on the desired resolution
        Msprev = M[τ-1,1:3];
        for s=1:Nsubsectors
            #generate pointcloud
            Ms = M[τ-1,1:3] + svec.*lineLength[τ].*s./Nsubsectors;  #nozzle position
            global Δt = (t[τ]-t[τ-1])./Nsubsectors;
            for y=Ms[3]-layerHeight[τ]:meshSize:Ms[3]
                #left and right sides
                if τ==2 || abs(lineFlow[τ-1])<=1e-9
                    append!(pointcloud, [Msprev[1],y]);
                    append!(normals, [-1,0]);
                    append!(boundaryNodes,size(pointcloud,2));
                    append!(depositionTime,t[τ-1]+(t[τ]-t[τ-1]).*s./Nsubsectors);
                end
                append!(pointcloud, [Ms[1],y]);
                append!(normals, [1,0]);
                append!(boundaryNodes,size(pointcloud,2));
                append!(depositionTime,t[τ-1]+(t[τ]-t[τ-1]).*s./Nsubsectors);
            end
            for x=Msprev[1]+meshSize:meshSize:Ms[1]-meshSize
                #bottom and top sides
                if layers[τ]==0
                    append!(pointcloud, [x,Ms[3]-layerHeight[τ]]);
                    append!(normals, [0,-1]);
                    append!(boundaryNodes,size(pointcloud,2));
                    append!(depositionTime,t[τ-1]+(t[τ]-t[τ-1]).*s./Nsubsectors);
                end
                append!(pointcloud, [x,Ms[3]]);
                append!(normals, [0,1]);
                append!(boundaryNodes,size(pointcloud,2));
                append!(depositionTime,t[τ-1]+(t[τ]-t[τ-1]).*s./Nsubsectors);
            end
            for x=Msprev[1]+meshSize:meshSize:Ms[1]-meshSize
                #internal nodes
                for y=Ms[3]-layerHeight[τ]+meshSize:meshSize:Ms[3]-meshSize
                    append!(pointcloud, [x,y]+(rand(Float64,2).-0.5).*meshSize/5);
                    append!(normals, [0,0]);
                    append!(internalNodes,size(pointcloud,2));
                    append!(depositionTime,t[τ-1]+(t[τ]-t[τ-1]).*s./Nsubsectors);
                end
            end
            Msprev = Ms;


            #neighbors search
            global Nprev = length(Nneighbors);  #number of nodes at the previous (sub)step
            global N = size(pointcloud,2);      #number of nodes at the current (sub)step
            global iter += 1;
            append!(T,T0*ones(N-Nprev));
            editedNodes = collect(Nprev+1:N);
            for i=length(Nneighbors)+1:N
                #initialize variables
                push!(neighbors,Int[]);
                push!(Nneighbors,0);
            end
            for i=Nprev+1:N
                #search the neighbors of the new nodes
                searchradius = minSearchRadius;
                while Nneighbors[i]<minNeighbors
                    #check every other node
                    for j=1:N
                        if j!=i && (pointcloud[1,j]-pointcloud[1,i])^2+(pointcloud[2,j]-pointcloud[2,i])^2<=searchradius^2
                            push!(neighbors[i],j);
                            #if the j-th node has been deposited previously
                            #add the node i to the j-th node neighbors list
                            if j<=Nprev && !(i in neighbors[j])
                                #println("Node ",i," added as new neighbor to the previous node ",j);
                                push!(neighbors[j],i);
                                #unique!(neighbors[j]);
                                #Nneighbors[j] = length(neighbors[j]);
                                Nneighbors[j] += 1;
                                push!(w[j],0);
                                push!(editedNodes,j);
                            end
                        end
                    end
                    unique!(neighbors[i]);
                    Nneighbors[i] = length(neighbors[i]);
                    searchradius += 0.5*minSearchRadius;
                end
            end

            #check boundary nodes
            for idx=1:lastindex(boundaryNodes)
                i = boundaryNodes[idx];
                for j in neighbors[i]
                    if dot(pointcloud[:,j]-pointcloud[:,i],normals[:,i])>1e-9
                        #the node i is no longer a boundary node
                        boundaryNodes[idx] = 0;
                        append!(internalNodes,i);
                        normals[:,i] = [0,0];
                    end
                end
            end
            deleteat!(boundaryNodes,findall(boundaryNodes.==0));

            #pointcloud plot
            #=if mod(s,3)==0
            #if layers[τ]==2
                figure();
                plot(pointcloud[1,boundaryNodes],pointcloud[2,boundaryNodes],"r.");
                plot(pointcloud[1,internalNodes],pointcloud[2,internalNodes],"k.");
                title("Pointcloud");
                axis("equal");
                display(gcf());
                return;
            end=#

            #neighbors distances and weights
            for i=Nprev+1:N
                #initialize variables
                push!(w,Array{Float64}(undef,Nneighbors[i]));
                push!(r2max,0);
            end
            for i in editedNodes
                r2 = zeros(Nneighbors[i]);
                for j=1:Nneighbors[i]
                    #P[i][:,j] = pointcloud[:,neighbors[i][j]]-pointcloud[:,i];
                    Prel = pointcloud[:,neighbors[i][j]]-pointcloud[:,i];
                    r2[j] = Prel'Prel;
                end
                r2max[i] = maximum(r2);
                for j=1:Nneighbors[i]
                    w[i][j] = exp(-6*r2[j]/r2max[i]);
                    #w[i][j] = 1.0;
                end
            end


            #thermal - boundary conditions
            for i=Nprev+1:N
                #initialize variables
                push!(g1,0);
                push!(g2,0);
                push!(g3,0);
            end
            for i in boundaryNodes
                if pointcloud[2,i]<=0+1e-6
                    #bottom surface
                    g1[i] = 1.0;
                    g2[i] = 0.0;
                    g3[i] = Tbed;
                else
                    #everywhere else
                    g1[i] = h;
                    g2[i] = k;
                    g3[i] = h*Tamb;
                end
            end
            #Fx = zeros(Float64,N);      #volumetric loads
            #Fy = zeros(Float64,N);


            #thermal - least square matrix inversion
            for i=Nprev+1:N
                #initialize variables
                push!(C,Matrix{Float64}(undef,6,Nneighbors[i]+2));
            end
            for i in editedNodes
                Prel = pointcloud[:,neighbors[i]].-pointcloud[:,i];
                xj = Prel[1,:];
                yj = Prel[2,:];
                if i in internalNodes
                    #internal node
                    V = zeros(Float64,1+Nneighbors[i],6);
                    for j=1:Nneighbors[i]
                        V[j,:] = [1, xj[j], yj[j], xj[j]^2, yj[j]^2, xj[j]*yj[j]];
                    end
                    V[1+Nneighbors[i],:] = [0, 0, 0, 2*k/(ρ*cp), 2*k/(ρ*cp), 0];
                    W = Diagonal(vcat(w[i],wpde));
                    (Q,R) = qr(W*V);
                    C[i] = inv(R)*transpose(Matrix(Q))*W;
                else
                    #boundary node
                    V = zeros(Float64,2+Nneighbors[i],6);
                    for j=1:Nneighbors[i]
                        V[j,:] = [1, xj[j], yj[j], xj[j]^2, yj[j]^2, xj[j]*yj[j]];
                    end
                    V[1+Nneighbors[i],:] = [0, 0, 0, 2*k/(ρ*cp), 2*k/(ρ*cp), 0];
                    V[2+Nneighbors[i],:] = [g1[i], g2[i]*normals[1,i], g2[i]*normals[2,i], 0, 0, 0];
                    W = Diagonal(vcat(w[i],wpde,wbc));
                    (Q,R) = qr(W*V);
                    C[i] = inv(R)*transpose(Matrix(Q))*W;
                end
            end


            #thermal - matrix assembly
            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
            Mt = sparse(rows,cols,vals,N,N);
        
        
            #thermal - linear system solution
            b = zeros(N);       #rhs vector
            for i=1:N
                b[i] = -C[i][1,1+Nneighbors[i]]*T[i]/Δt;
            end
            for i in boundaryNodes
                b[i] += C[i][1,2+Nneighbors[i]]*g3[i];
            end
            Tnew = Mt\b;
            dTdt = (Tnew-T)./Δt;
            global T = Tnew;

            
            #heat fluxes
            global dTdx = zeros(N);
            global dTdy = zeros(N);
            for i=1:N
                for j=1:Nneighbors[i]
                    global dTdx[i] += C[i][2,j]*T[neighbors[i][j]];        #C[i][j]*u[j]
                    global dTdy[i] += C[i][3,j]*T[neighbors[i][j]];
                end
                global dTdx[i] += C[i][2,1+Nneighbors[i]]*dTdt[i];     #C[i][1+N]*gPDE
                global dTdy[i] += C[i][3,1+Nneighbors[i]]*dTdt[i];
            end
            for i in boundaryNodes
                global dTdx[i] += C[i][2,2+Nneighbors[i]]*g3[i];       #C[i][2+N]*gBC
                global dTdy[i] += C[i][3,2+Nneighbors[i]]*g3[i];
            end
            global qx = -k*dTdx;
            global qy = -k*dTdy;

            #mechanical - boundary conditions
            for i=Nprev+1:N
                #initialize variables
                push!(g1u,0);
                push!(g2u,0);
                push!(g3u,0);
                push!(g1v,0);
                push!(g2v,0);
                push!(g3v,0);
            end
            for i in boundaryNodes
                if pointcloud[2,i]<=0+1e-6
                    #bottom surface
                    g1u[i] = 1.0;
                    g2u[i] = 0.0;
                    g3u[i] = 0.0;
                    g1v[i] = 1.0;
                    g2v[i] = 0.0;
                    g3v[i] = 0.0;
                else
                    #everywhere else
                    g1u[i] = 0.0;
                    g2u[i] = 1.0;
                    g3u[i] = 0.0;
                    g1v[i] = 0.0;
                    g2v[i] = 1.0;
                    g3v[i] = 0.0;
                end
            end
            Fx = zeros(Float64,N);
            Fy = zeros(Float64,N);

            #mechanical - least square matrix inversion
            for i=Nprev+1:N
                #initialize variables
                push!(Cm,Matrix{Float64}(undef,6,Nneighbors[i]+2));
            end
            for i in editedNodes
                Prel = pointcloud[:,neighbors[i]].-pointcloud[:,i];
                xj = Prel[1,:];
                yj = Prel[2,:];
                if i in internalNodes
                    #internal node
                    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*1/(1-ν^2), 2*0.5*1/(1+ν), 0,  0, 0, 0, 0, 0, 1*ν/(1-ν^2)+0.5*1/(1+ν)];
                    V[2+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 1*ν/(1-ν^2)+0.5*1/(1+ν),  0, 0, 0, 2*0.5*1/(1+ν), 2*1/(1-ν^2), 0];
                    W = Diagonal(vcat(w[i],w[i],wpde,wpde));
                    VF = svd(W*V);
                    Cm[i] = transpose(VF.Vt)*inv(Diagonal(VF.S))*transpose(VF.U)*W;
                else
                    #boundary node
                    nx = normals[1,i];
                    ny = normals[2,i];
                    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*1/(1-ν^2), 2*0.5*1/(1+ν), 0,  0, 0, 0, 0, 0, 1*ν/(1-ν^2)+0.5*1/(1+ν)];
                    V[2+2*Nneighbors[i],:] = [0, 0, 0, 0, 0, 1*ν/(1-ν^2)+0.5*1/(1+ν),  0, 0, 0, 2*0.5*1/(1+ν), 2*1/(1-ν^2), 0];
                    V[3+2*Nneighbors[i],:] = [g1u[i], g2u[i]*nx*1/(1-ν^2), g2u[i]*ny*0.5*1/(1+ν), 0, 0, 0,  0, g2u[i]*ny*0.5*1/(1+ν), g2u[i]*nx*1*ν/(1-ν^2), 0, 0, 0];
                    V[4+2*Nneighbors[i],:] = [0, g2v[i]*ny*1*ν/(1-ν^2), g2v[i]*nx*0.5*1/(1+ν), 0, 0, 0,  g1v[i], g2v[i]*nx*0.5*1/(1+ν), g2v[i]*ny*1/(1-ν^2), 0, 0, 0];
                    W = Diagonal(vcat(w[i],w[i],wpde,wpde,wbc,wbc));
                    VF = svd(W*V);
                    Cm[i] = transpose(VF.Vt)*inv(Diagonal(VF.S))*transpose(VF.U)*W;
                end
            end


            #mechanical - matrix assembly
            rowsM = Int[];
            colsM = Int[];
            valsM = Float64[];
            for i=1:N
                #u equation
                push!(rowsM, i);
                push!(colsM, i);
                push!(valsM, 1);
                for j=1:Nneighbors[i]
                    push!(rowsM, i);
                    push!(colsM, neighbors[i][j]);
                    push!(valsM, -Cm[i][1,j]);
                    push!(rowsM, i);
                    push!(colsM, N+neighbors[i][j]);
                    push!(valsM, -Cm[i][1,j+Nneighbors[i]]);
                end

                #v equation
                push!(rowsM, N+i);
                push!(colsM, N+i);
                push!(valsM, 1);
                for j=1:Nneighbors[i]
                    push!(rowsM, N+i);
                    push!(colsM, neighbors[i][j]);
                    push!(valsM, -Cm[i][7,j]);
                    push!(rowsM, N+i);
                    push!(colsM, N+neighbors[i][j]);
                    push!(valsM, -Cm[i][7,j+Nneighbors[i]]);
                end
            end
            Mm = sparse(rowsM,colsM,valsM,2*N,2*N);


            #mechanical - linear system solution
            bm = zeros(2*N);       #rhs vector
            for i=1:N
                bm[i] = 0;
                bm[N+i] = 0;
                bm[i] += Cm[i][1,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][1,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
                bm[N+i] += Cm[i][7,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][7,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
            end
            for i in boundaryNodes
                nx = normals[1,i];
                ny = normals[2,i];
                bm[i] += Cm[i][1,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][1,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));
                bm[N+i] += Cm[i][7,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][7,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));
            end
            sol = Mm\bm;
            global u = sol[1:N];
            global v = sol[N+1:2*N];

            
            #stresses
            global dudx = zeros(N);
            global dudy = zeros(N);
            global dvdx = zeros(N);
            global dvdy = zeros(N);
            for i=1:N
                for j=1:Nneighbors[i]
                    global dudx[i] += Cm[i][2,j]*u[neighbors[i][j]];        #C[i][j]*u[j]
                    global dudy[i] += Cm[i][3,j]*u[neighbors[i][j]];
                    global dvdx[i] += Cm[i][8,j]*v[neighbors[i][j]];        #C[i][j]*v[j]
                    global dvdy[i] += Cm[i][9,j]*v[neighbors[i][j]];
                end
                global dudx[i] += Cm[i][2,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][2,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));     #C[i][1+N]*gPDE
                global dudy[i] += Cm[i][3,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][3,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
                global dvdx[i] += Cm[i][8,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][8,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));     #C[i][1+N]*gPDE
                global dvdy[i] += Cm[i][9,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][9,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
            end
            for i in boundaryNodes
                nx = normals[1,i];
                ny = normals[2,i];
                global dudx[i] += Cm[i][2,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][2,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));       #C[i][2+N]*gBC
                global dudy[i] += Cm[i][3,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][3,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));
            end
            global σxx = @. (E/(1-ν^2))*(1*dudx + ν*dvdy) - E*α*(T-Tref)/(1-ν);
            global σyy = @. (E/(1-ν^2))*(ν*dudx + 1*dvdy) - E*α*(T-Tref)/(1-ν);
            global τxy = @. (0.5*E/(1+ν))*(dudy + dvdx);
            global σv = @. sqrt(σxx^2 + σyy^2 - σxx*σyy + 3*τxy^2);



            #deposition time plot
            #=if mod(s,10)==0
                figure();
                scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=depositionTime,cmap="inferno");
                colorbar();
                title("Deposition time (s)");
                axis("equal");
                display(gcf());
                return;
            end=#

            #connectivity plot
            #=if mod(s,2)==0
                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());
                return;
            end=#


            #transient plot - mechanical
            #=if mod(s,10)==0
                figure();
                scatter(pointcloud[1,:],pointcloud[2,:],c=sqrt.(u.^2+v.^2),cmap="Oranges");
                title("Displacement - t = "*string(t[τ-1]+Δt));
                colorbar();
                axis("equal");
                display(gcf());
                return;
            end=#
            
            if savePlots
                #temperature plot
                figure();
                plt = scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=T,cmap="inferno");
                title("Temperature - t = "*string(round(depositionTime[end],digits=5)));
                plt.set_clim(Tbed-10,T0+10);
                colorbar(plt);
                xlabel("X coordinate (mm)");
                ylabel("Z coordinate (mm)");
                axis("equal");
                savefig("vid1/"*string(lpad(iter,3,"0"))*".jpg");
                #display(gcf());

                #stresses plot
                figure();
                pltv = scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=σv,cmap="jet");
                title("Von Mises stress - t = "*string(round(depositionTime[end],digits=5)));
                pltv.set_clim(0,1e9);
                colorbar(pltv);
                xlabel("X coordinate (mm)");
                ylabel("Z coordinate (mm)");
                axis("equal");
                savefig("vid2/"*string(lpad(iter,3,"0"))*".jpg");
                #display(gcf());
            end
        end
    else
        #rapid movement without extrusion
    end
end

#cooling after print
Δt = 10;
tcooling = 60;
ti = t[end];
tf = ti+tcooling;
while ti<tf
    global ti += Δt;

    #thermal - matrix assembly
    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
    Mt = sparse(rows,cols,vals,N,N);


    #thermal - linear system solution
    b = zeros(N);       #rhs vector
    for i=1:N
        b[i] = -C[i][1,1+Nneighbors[i]]*T[i]/Δt;
    end
    for i in boundaryNodes
        b[i] += C[i][1,2+Nneighbors[i]]*g3[i];
    end
    Tnew = Mt\b;
    dTdt = (Tnew-T)./Δt;
    global T = Tnew;


    #heat fluxes
    global dTdx = zeros(N);
    global dTdy = zeros(N);
    for i=1:N
        for j=1:Nneighbors[i]
            global dTdx[i] += C[i][2,j]*T[neighbors[i][j]];        #C[i][j]*u[j]
            global dTdy[i] += C[i][3,j]*T[neighbors[i][j]];
        end
        global dTdx[i] += C[i][2,1+Nneighbors[i]]*dTdt[i];     #C[i][1+N]*gPDE
        global dTdy[i] += C[i][3,1+Nneighbors[i]]*dTdt[i];
    end
    for i in boundaryNodes
        global dTdx[i] += C[i][2,2+Nneighbors[i]]*g3[i];       #C[i][2+N]*gBC
        global dTdy[i] += C[i][3,2+Nneighbors[i]]*g3[i];
    end
    global qx = -k*dTdx;
    global qy = -k*dTdy;


    #mechanical - matrix assembly
    Fx = zeros(Float64,N);
    Fy = zeros(Float64,N);
    rowsM = Int[];
    colsM = Int[];
    valsM = Float64[];
    for i=1:N
        #u equation
        push!(rowsM, i);
        push!(colsM, i);
        push!(valsM, 1);
        for j=1:Nneighbors[i]
            push!(rowsM, i);
            push!(colsM, neighbors[i][j]);
            push!(valsM, -Cm[i][1,j]);
            push!(rowsM, i);
            push!(colsM, N+neighbors[i][j]);
            push!(valsM, -Cm[i][1,j+Nneighbors[i]]);
        end

        #v equation
        push!(rowsM, N+i);
        push!(colsM, N+i);
        push!(valsM, 1);
        for j=1:Nneighbors[i]
            push!(rowsM, N+i);
            push!(colsM, neighbors[i][j]);
            push!(valsM, -Cm[i][7,j]);
            push!(rowsM, N+i);
            push!(colsM, N+neighbors[i][j]);
            push!(valsM, -Cm[i][7,j+Nneighbors[i]]);
        end
    end
    Mm = sparse(rowsM,colsM,valsM,2*N,2*N);

    
    #mechanical - linear system solution
    bm = zeros(2*N);       #rhs vector
    for i=1:N
        bm[i] = 0;
        bm[N+i] = 0;
        bm[i] += Cm[i][1,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][1,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
        bm[N+i] += Cm[i][7,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][7,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
    end
    for i in boundaryNodes
        nx = normals[1,i];
        ny = normals[2,i];
        bm[i] += Cm[i][1,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][1,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));
        bm[N+i] += Cm[i][7,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][7,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));
    end
    sol = Mm\bm;
    global u = sol[1:N];
    global v = sol[N+1:2*N];
    
    
    #stresses
    global dudx = zeros(N);
    global dudy = zeros(N);
    global dvdx = zeros(N);
    global dvdy = zeros(N);
    for i=1:N
        for j=1:Nneighbors[i]
            global dudx[i] += Cm[i][2,j]*u[neighbors[i][j]];        #C[i][j]*u[j]
            global dudy[i] += Cm[i][3,j]*u[neighbors[i][j]];
            global dvdx[i] += Cm[i][8,j]*v[neighbors[i][j]];        #C[i][j]*v[j]
            global dvdy[i] += Cm[i][9,j]*v[neighbors[i][j]];
        end
        global dudx[i] += Cm[i][2,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][2,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));     #C[i][1+N]*gPDE
        global dudy[i] += Cm[i][3,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][3,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
        global dvdx[i] += Cm[i][8,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][8,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));     #C[i][1+N]*gPDE
        global dvdy[i] += Cm[i][9,1+2*Nneighbors[i]]*(-Fx[i]/E + 1*α*dTdx[i]/(1-ν)) + Cm[i][9,2+2*Nneighbors[i]]*(-Fy[i]/E + 1*α*dTdy[i]/(1-ν));
    end
    for i in boundaryNodes
        nx = normals[1,i];
        ny = normals[2,i];
        global dudx[i] += Cm[i][2,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][2,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));       #C[i][2+N]*gBC
        global dudy[i] += Cm[i][3,3+2*Nneighbors[i]]*(g3u[i]/E + g2u[i]*nx*1*α*(T[i]-Tref)/(1-ν)) + Cm[i][3,4+2*Nneighbors[i]]*(g3v[i]/E + g2v[i]*ny*1*α*(T[i]-Tref)/(1-ν));
    end
    global σxx = @. (E/(1-ν^2))*(1*dudx + ν*dvdy) - E*α*(T-Tref)/(1-ν);
    global σyy = @. (E/(1-ν^2))*(ν*dudx + 1*dvdy) - E*α*(T-Tref)/(1-ν);
    global τxy = @. (0.5*E/(1+ν))*(dudy + dvdx);
    global σv = @. sqrt(σxx^2 + σyy^2 - σxx*σyy + 3*τxy^2);


    
    if savePlots
        #temperature plot
        figure();
        plt = scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=T,cmap="inferno");
        title("Temperature - t = "*string(round(t[τ],digits=5)));
        plt.set_clim(Tbed-10,T0+10);
        colorbar(plt);
        xlabel("X coordinate (mm)");
        ylabel("Z coordinate (mm)");
        axis("equal");
        savefig("vid1/"*string(lpad(iter,3,"0"))*".jpg");
        #display(gcf());

        #stresses plot
        figure();
        pltv = scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=σv,cmap="jet");
        title("Von Mises stress - t = "*string(round(t[τ],digits=5)));
        pltv.set_clim(0,1e9);
        colorbar(pltv);
        xlabel("X coordinate (mm)");
        ylabel("Z coordinate (mm)");
        axis("equal");
        savefig("vid2/"*string(lpad(iter,3,"0"))*".jpg");
        #display(gcf());
    end
end


println("Simulation completed in ",round(time()-time0,digits=2)," s");
println("Maximum stress: ",maximum(σv)," at index ",findfirst(σv.==maximum(σv)));

#deposition time plot
figure();
scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=depositionTime,cmap="inferno");
colorbar();
title("Deposition time (s)");
xlabel("X coordinate (mm)");
ylabel("Z coordinate (mm)");
axis("equal");
display(gcf());

#temperature plot
figure();
scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=T,cmap="inferno");
title("Numerical solution - temperature");
xlabel("X coordinate (mm)");
ylabel("Z coordinate (mm)");
colorbar();
axis("equal");
display(gcf());

#displacement plot
figure();
scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=sqrt.(u.^2+v.^2),cmap="Oranges");
title("Numerical solution - displacement");
xlabel("X coordinate (mm)");
ylabel("Z coordinate (mm)");
colorbar();
axis("equal");
display(gcf());

#von Mises stress plot
figure();
pltv = scatter(pointcloud[1,:]*1000,pointcloud[2,:]*1000,c=σv,cmap="jet");
title("Numerical solution - von Mises stress");
pltv.set_clim(0,1e9);
colorbar(pltv);
xlabel("X coordinate (mm)");
ylabel("Z coordinate (mm)");
axis("equal");
display(gcf());
