mesh.jl

export init_grid, save

using WriteVTK

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

Base.show(io::IO, ::MIME"text/plain", f::Mesh) =
    print("$(nameof(typeof(f))), $(size(f.cells, 2)) cells")

function init_grid(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 = [x[i] for i in 1:2, x in vec(coords)]
    cells = Array{Int, 2}(undef, 3, 2 * m * n)

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

    return Mesh(vertices, cells)
end

init_grid(img::Array{<:Any, 2}, type=:vertex) =
    type == :vertex ?
	init_grid(size(img, 1) - 1, size(img, 2) - 1, (1.0, 1.0), size(img)) :
	init_grid(size(img, 1), size(img, 2), (0.5, 0.5), size(img) .- (0.5, 0.5))

#function cell_contains(mesh, cell, v)
#    geo = mesh.vertices[:, mesh.cells[:, cell]]
#    J = jacobian(x -> geo * [1 - x[1] - x[2], x[1], x[2]], [0., 0.])
#    λ = J \ (v - geo[:, 1])
#    return all(λ .>= 0) && sum(λ) <= 1
#end

function vtk_mesh(filename, mesh::Mesh)
    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, 3*(i-1)+1:3*(i-1)+3)
	for i in axes(mesh.cells, 2)]
    vertices = reshape(mesh.vertices[:, mesh.cells], 2, :)

    return vtk_grid(filename, vertices, cells)
end

"convenience function for saving to vtk"
function save(filename, mesh::Mesh, fs...)
    vtk = vtk_mesh(filename, mesh)
    for f in fs
	f.space.mesh == mesh ||
	    throw(ArgumentError("meshes do not match"))
        vtk_append!(vtk, f)
    end
    vtk_save(vtk)
end