#----------------------------------------------------------------------
#  GFDM method test - 3d heat equation, Mixed BC
#  Author: Andrea Pavan
#  Date: 22/12/2022
#  License: GPLv3-or-later
#----------------------------------------------------------------------
using ElasticArrays;
using LinearAlgebra;
using SparseArrays;
using Printf;
using PyPlot;
include("utils.jl");


#problem definition
l1 = 5.0;       #domain x size
l2 = 2.5;       #domain y size
l3 = 2.0;       #domain z size
uL = 400;       #left border temperature
uR = 300;       #right border temperature

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


#pointcloud generation - boundary nodes
time1 = time();
pointcloud = ElasticArray{Float64}(undef,3,0);      #3xN matrix containing the coordinates [X;Y;Z] of each node
boundaryNodes = Vector{Int};        #indices of the boundary nodes
normals = ElasticArray{Float64}(undef,3,0);     #3xN matrix containing the components [nx;ny;nz] of the normal of each boundary node
(section,sectionnormals) = defaultCrossSection(l2, l3, surfaceMeshSize);
for x in 0:surfaceMeshSize:l1
    append!(pointcloud, vcat(zeros(Float64,1,size(section,2)).+x,section));
    append!(normals, vcat(zeros(Float64,1,size(sectionnormals,2)),sectionnormals));
end
for y in -l2/2+surfaceMeshSize:surfaceMeshSize:l2/2-surfaceMeshSize
    for z in -l3+surfaceMeshSize:surfaceMeshSize:0-surfaceMeshSize
        if abs(y)<(l2-l3)/2 || (abs(y)-(l2-l3)/2)^2+(z+l3/2)^2<(l3/2)^2
            append!(pointcloud, [0,y,z]);
            append!(normals, [-1,0,0]);
            append!(pointcloud, [l1,y,z]);
            append!(normals, [1,0,0]);
        end
    end
end
boundaryNodes = collect(range(1,size(pointcloud,2)));


#pointcloud generation - internal nodes (linspaced)
#=surfaceMargin = surfaceMeshSize/2;
for y in -l2/2+surfaceMargin:meshSize:l2/2-surfaceMargin
    for z in -l3+surfaceMargin:meshSize:0-surfaceMargin
        if abs(y)<(l2-l3)/2 || (abs(y)-(l2-l3)/2)^2+(z+l3/2)^2<(l3/2)^2
            for x in 0+meshSize:meshSize:l1-meshSize
                append!(pointcloud, [x,y,z]+(rand(Float64,3).-0.5).*meshSize/5);
            end
        end
    end
end=#


#pointcloud generation - internal nodes (octree)
#=(octree,octreeSize,octreeCenter,octreePoints,octreeNpoints) = buildOctree(pointcloud);
octreeLeaves = findall(octreeNpoints.>=0);
for i in octreeLeaves
    if octreeNpoints[i] != 1
        if (abs(octreeCenter[2,i])<(l2-l3)/2 && octreeCenter[3,i]<=0) || (abs(octreeCenter[2,i])-(l2-l3)/2)^2+(octreeCenter[3,i]+l3/2)^2<(l3/2)^2
            append!(pointcloud, octreeCenter[:,i]);
        end
    end
end=#


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


#generate cartesian subgrid
N = size(pointcloud,2);     #number of points
Pmax = maximum(pointcloud,dims=2);
Pmin = minimum(pointcloud,dims=2);
boundingBoxSize = Pmax-Pmin;
Ncellx = max(1,Int(round(boundingBoxSize[1]/(3*meshSize))));
Ncelly = max(1,Int(round(boundingBoxSize[2]/(3*meshSize))));
Ncellz = max(1,Int(round(boundingBoxSize[3]/(3*meshSize))));
cellsizex = boundingBoxSize[1]/Ncellx;
cellsizey = boundingBoxSize[2]/Ncelly;
cellsizez = boundingBoxSize[3]/Ncellz;
cell = Array{Vector{Int},3}(undef,Ncellx,Ncelly,Ncellz);
cellIdx = Matrix{Int}(undef,3,N);

function repopulateCells()
    for i=1:Ncellx
        for j=1:Ncelly
            for k=1:Ncellz
                global cell[i,j,k] = Vector{Int}(undef,0);
            end
        end
    end
    global cellIdx = Matrix{Int}(undef,3,N);
    for i=1:N
        idx1 = max(1,min(Ncellx,1+Int(floor((pointcloud[1,i]-Pmin[1])/cellsizex))));
        idx2 = max(1,min(Ncelly,1+Int(floor((pointcloud[2,i]-Pmin[2])/cellsizey))));
        idx3 = max(1,min(Ncellz,1+Int(floor((pointcloud[3,i]-Pmin[3])/cellsizez))));
        push!(cell[idx1,idx2,idx3], i);
        global cellIdx[:,i] = [idx1,idx2,idx3];
    end
end

#pointcloud generation - internal nodes (flowmesher)
v = zeros(Float64,3,N);
F = zeros(Float64,3,N);
maxV = 1.0;
dt = 0.1;
h = 3*meshSize;
iter = 1;
while maxV>5e-3
    repopulateCells();
    for i=1+length(boundaryNodes):N
        global F[:,i] = [0,0,0];
        for j=1:N
            rd2 = transpose(pointcloud[:,j]-pointcloud[:,i])*(pointcloud[:,j]-pointcloud[:,i]);
            if i!=j && rd2<h^2
                global F[:,i] += exp(-6.5*rd2/(h^2))*(pointcloud[:,i]-pointcloud[:,j])/sqrt(rd2);
            end
        end
        #=for j in cell[cellIdx[1,i],cellIdx[2,i],cellIdx[3,i]];
            rd2 = transpose(pointcloud[:,j]-pointcloud[:,i])*(pointcloud[:,j]-pointcloud[:,i]);
            if i!=j && rd2<h^2
                global F[:,i] += exp(-6.5*rd2/(h^2))*(pointcloud[:,i]-pointcloud[:,j])/sqrt(rd2);
            end
        end=#
        #=cellstosearch = [cellIdx[:,i],
                        [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                        [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                        [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                        [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                        [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], cellIdx[3,i]],
                        [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                        [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                        [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                        [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)],
                        
                        [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                        [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                        [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                        [cellIdx[1,i], cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                        [cellIdx[1,i], cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                        [cellIdx[1,i], max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                        [cellIdx[1,i], max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                        [cellIdx[1,i], max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)],

                        [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                        [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                        [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                        [max(1,cellIdx[1,i]-1), cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                        [max(1,cellIdx[1,i]-1), cellIdx[2,i], cellIdx[3,i]],
                        [max(1,cellIdx[1,i]-1), cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                        [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                        [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                        [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)]];
        unique!(cellstosearch);
        for cindex in cellstosearch
            #find the neighbors on cell[cindex]
            for j in cell[cindex[1],cindex[2],cindex[3]]
                rd2 = transpose(pointcloud[:,j]-pointcloud[:,i])*(pointcloud[:,j]-pointcloud[:,i]);
                if i!=j && rd2<h^2
                    global F[:,i] += exp(-6.5*rd2/(h^2))*(pointcloud[:,i]-pointcloud[:,j])/sqrt(rd2);
                end
            end
        end=#
    end
    global v += dt*F;
    for i=1+length(boundaryNodes):N
        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 pointcloud += dt.*v;
    for i=1+length(boundaryNodes):N
        #if abs(y)<(l2-l3)/2 || (abs(y)-(l2-l3)/2)^2+(z+l3/2)^2<(l3/2)^2
        if pointcloud[3,i]>0-surfaceMargin && abs(pointcloud[2,i])<(l2-l3)/2
            pointcloud[3,i] -= 1.5*(pointcloud[3,i]+surfaceMargin);
            v[:,i] = [0,0,0];
        end
        if pointcloud[3,i]<-l2+surfaceMargin && abs(pointcloud[2,i])<(l2-l3)/2
            pointcloud[3,i] -= 1.5*(pointcloud[3,i]-l2+surfaceMargin);
            v[:,i] = [0,0,0];
        end
        if pointcloud[1,i]>l1-surfaceMargin && abs(pointcloud[2,i])<(l2-l3)/2
            pointcloud[1,i] -= 1.5*(pointcloud[1,i]-l1+surfaceMargin);
            v[:,i] = [0,0,0];
        end
        if pointcloud[1,i]<0+surfaceMargin && abs(pointcloud[2,i])<(l2-l3)/2
            pointcloud[1,i] -= 1.5*(pointcloud[1,i]+0-surfaceMargin);
            v[:,i] = [0,0,0];
        end
        if abs(pointcloud[2,i])>=(l2-l3)/2 && (abs(pointcloud[2,i])-(l2-l3)/2)^2+(pointcloud[3,i]+l3/2)^2>(l3/2-surfaceMargin)^2
            relpos = pointcloud[2:3,i]-[sign(pointcloud[2,i])*(l2-l3)/2, -l3/2];
            pointcloud[2:3,i] -= relpos./2;
            v[:,i] = [0,0,0];
        end
        #=if pointcloud[1,i]^2+pointcloud[2,i]^2<(a+surfaceMargin)^2
            ang = atan(pointcloud[2,i],pointcloud[1,i]);
            rad = sqrt(pointcloud[:,i]'pointcloud[:,i]);
            pointcloud[:,i] += 2*((a+surfaceMargin).*[cos(ang),sin(ang)]-pointcloud[:,i]);
            v[:,i] = [0,0];
        end=#
    end
    global iter += 1;
    if mod(iter,10)==0
        figure();
        plot(pointcloud[1,1+length(boundaryNodes):N],pointcloud[2,1+length(boundaryNodes):N],"k.");
        #plot(pointcloud[2,1+length(boundaryNodes):N],pointcloud[3,1+length(boundaryNodes):N],"k.");
        title("Pointcloud plot");
        axis("equal");
        display(gcf());
        sleep(0.1);
    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));

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