runtests.jl

import Random
using Test

using SemiSmoothNewton

Random.seed!(0)

function test_mesh(mesh)
    # assemble edge -> cells map
    edgemap = Dict{NTuple{2, Int}, Vector{Int}}()
    for cell in cells(mesh)
	vs = sort(SVector(vertices(mesh, cell)))

	e1 = (vs[1], vs[2])
	e2 = (vs[1], vs[3])
	e3 = (vs[2], vs[3])

	edgemap[e1] = push!(get!(edgemap, e1, []), cell)
	edgemap[e2] = push!(get!(edgemap, e2, []), cell)
	edgemap[e3] = push!(get!(edgemap, e3, []), cell)
    end
    for (edge, cells) in edgemap
	@assert(1 <= length(cells) <= 2)
    end
    # are cells positively oriented?
    for cell in cells(mesh)
	delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
	@assert(det(delmap) > 0)
    end
end

@testset "imaging" begin
    img = rand(3, 3)
    mesh = init_grid(img)

    @test ndims_domain(mesh) == 2
    @test ndims_space(mesh) == 2

    # dual mesh: vertices are pixel centers
    @test nvertices(mesh) == 3^2
    @test ncells(mesh) == 2 * 2^2 # TODO: only simplex mesh
    v0 = minimum(mesh.vertices, dims = 2) |> vec
    v1 = maximum(mesh.vertices, dims = 2) |> vec
    @test v0 == [1., 1.]
    @test v1 == [3., 3.]

    u_space = FeSpace(mesh, P1(), (1,))
    u = FeFunction(u_space)

    interpolate!(u, img)
    img_sampled = sample(u)

    @test img == img_sampled
end

@testset "mesh refinement" begin
    # simple hashing
    val(f) = integrate(f.space.mesh, (x; f) -> sum(f .* f); f)

    mesh = init_grid(zeros(5, 5))

    f = FeFunction(FeSpace(mesh, P1(), (2, 3)))
    f.data .= rand(size(f.data)...)

    for i = 1:2
        # refine roughly half of all cells
        cells_ref = Set(findall(_ -> rand(Bool), cells(mesh)))

        mesh_new, (f_new,) = refine(mesh, cells_ref; f)

        @test isapprox(val(f), val(f_new))

        (mesh, f) = (mesh_new, f_new)
    end
end