mesh.jl

export init_grid, init_hgrid, save, refine!, refine, cells, vertices,
    ndims_domain, ndims_space, nvertices, nvertices_cell, ncells

using LinearAlgebra: norm

using WriteVTK
using StaticArrays: SVector

struct HMesh
    vertices::Vector{SVector{2, Float64}}
    cells::Vector{NTuple{3, Int}}
    levels::Vector{Int}
end

Base.show(io::IO, ::MIME"text/plain", x::HMesh) =
    print("$(nameof(typeof(x))), $(ncells(x)) cells")

ndims_domain(::HMesh) = 2
ndims_space(::HMesh) = 2
# TODO: move to simplex mesh interface
nvertices_cell(::HMesh) = 3
# TODO: do we even these two?
nvertices(x::HMesh) = length(x.vertices)
ncells(x::HMesh) = length(x.cells)

cells(mesh::HMesh) = setdiff(axes(mesh.cells, 1), findall(>(0), mesh.levels))
vertices(mesh::HMesh, cell) = mesh.cells[cell]

function init_hgrid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
    r1 = LinRange(v0[1], v1[1], m + 1)
    r2 = LinRange(v0[2], v1[2], n + 1)
    coords = collect(Iterators.product(r1, r2))

    vertices = SVector{2, Float64}.(vec(coords))
    cells = Vector{NTuple{3, Int}}()

    vidx = LinearIndices((m + 1, n + 1))
    e1 = CartesianIndex(1, 0)
    e2 = CartesianIndex(0, 1)
    k = 0
    for I in CartesianIndices((m, n))
	push!(cells, (vidx[I], vidx[I + e1], vidx[I + e1 + e2]))
	push!(cells, (vidx[I], vidx[I + e1 + e2], vidx[I + e2]))
    end

    return HMesh(vertices, cells, zeros(Int, axes(cells)))
end

function init_hgrid(img::Array{<:Any, 2}; type=:vertex)
    s = (size(img, 2), size(img, 1))
    type == :vertex ?
        init_hgrid((s .- 1)..., (1.0, 1.0), s) :
	init_hgrid(s..., (0.5, 0.5), s .- (0.5, 0.5))
end

function sub_mesh(hmesh::HMesh, subcells = cells(hmesh))
    cells = collect(reshape(reinterpret(Int, hmesh.cells[subcells]), 3, :))
    # TODO: restrict vertices too
    vertices = collect(reshape(reinterpret(Float64, hmesh.vertices), 2, :))
    return Mesh(vertices, cells)
end


# 2d, simplex grid
struct Mesh
    vertices::Array{Float64, 2}
    cells::Array{Int, 2}
end

# needs to be after Mesh definition