diff --git a/src/SemiSmoothNewton.jl b/src/SemiSmoothNewton.jl
index 7a73355b1145a1bd23853f5e5226a1e1d49c86d0..d52e24d12cc38c3b55f41a4a44d264f06451e82e 100644
--- a/src/SemiSmoothNewton.jl
+++ b/src/SemiSmoothNewton.jl
@@ -5,56 +5,4 @@ include("function.jl")
include("operator.jl")
include("image.jl")
-function clip_line(r0, r1, l0, l1)
- d_l = -(l1[1] - l0[1])
- d_r = +(l1[1] - l0[1])
- d_b = -(l1[2] - l0[2])
- d_t = +(l1[2] - l0[2])
-
- t_l = -(r0[1] - l0[1]) / d_l
- t_r = +(r1[1] - l0[1]) / d_r
- t_b = -(r0[2] - l0[2]) / d_b
- t_t = +(r1[2] - l0[2]) / d_t
-
- t0 = max(
- d_l <= 0 ? t_l : 0.,
- d_r <= 0 ? t_r : 0.,
- d_b <= 0 ? t_b : 0.,
- d_t <= 0 ? t_t : 0.)
-
- t1 = min(
- d_l >= 0 ? t_l : 1.,
- d_r >= 0 ? t_r : 1.,
- d_b >= 0 ? t_b : 1.,
- d_t >= 0 ? t_t : 1.)
-
- # tolerance on non-intersection
- t1 < t0 + eps() && return false
-
- c0 = l0 + t0 * (l1 - l0)
- c1 = l0 + t1 * (l1 - l0)
-
- return c0, c1
-end
-
-"project array data onto a grid function"
-function project_array!(u, img::Array)
- vertices = u.mesh.vertices
- mesh_box = (
- reduce((a, b) -> min.(a, b), vertices),
- reduce((a, b) -> max.(a, b), vertices))
-
-
- for cell in u.mesh.cells
- cell_box = (
- reduce((a, b) -> min.(a, b), cell.vertices),
- reduce((a, b) -> max.(a, b), cell.vertices))
-
- u[cell]
-
- end
-end
-
-
-
end # module
diff --git a/src/image.jl b/src/image.jl
index 10e3740e74bc7a52896bf6bf2b43800a39fdaa5b..b770b7fe3275d549b686b2cfef8288133ea0f41b 100644
--- a/src/image.jl
+++ b/src/image.jl
@@ -15,9 +15,6 @@ function peak_signal_to_noise_ratio(img, img_ref)
return 10 * log10(imax / mse)
end
-# FIXME: unused?
-from_img(img) = permutedims(reverse(arr; dims = 1))
-
eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]
# uses native array indexing, i.e. 1:n
diff --git a/src/mesh.jl b/src/mesh.jl
index fc82db5c8a3d988845e83078548aaa0bb6dd4972..a9d4eefc9c0f810360e8a4c18541f50eec3f174b 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -97,36 +97,6 @@ cells(mesh::Mesh) = axes(mesh.cells, 2)
vertices(mesh::Mesh) = axes(mesh.vertices, 2)
vertices(mesh::Mesh, cell) = ntuple(i -> mesh.cells[i, cell], nvertices_cell(mesh))
-#function Base.getproperty(obj::Mesh, sym::Symbol)
-# if sym === :vertices
-# return reinterpret(SVector{ndims_space(obj), Float64}, vec(obj.vertices))
-# else
-# return getfield(obj, sym)
-# end
-#end
-
-# old grid initialization
-# produces regular grid but not suitable for newest vertex bisection
-function init_grid_old(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
-
# new grid initialization
# produces regular grid suitable for newest vertex bisection
function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
@@ -355,13 +325,6 @@ function geo_contains(A, v)
return all(λ .>= 0) && sum(λ) <= 1
end
-#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)]
@@ -424,27 +387,3 @@ function elmap(mesh, cell)
@inbounds view(mesh.vertices, :, view(mesh.cells, :, cell)))
return x -> A * SA[1 - x[1] - x[2], x[1], x[2]]
end
-
-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
diff --git a/test/runtests.jl b/test/runtests.jl
index 3f7ca2d746492329c9c825a960a4fe977e6d7526..d3c8789b5a7dd84fc2dc9111417370292cd14909 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,30 @@ 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)