cleanup

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

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

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

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)