diff --git a/src/SemiSmoothNewton.jl b/src/SemiSmoothNewton.jl
index 47b41372fbecc66049c15a2184f0bfc212fa6425..7a73355b1145a1bd23853f5e5226a1e1d49c86d0 100644
--- a/src/SemiSmoothNewton.jl
+++ b/src/SemiSmoothNewton.jl
@@ -1,9 +1,9 @@
module SemiSmoothNewton
-include("image.jl")
include("mesh.jl")
include("function.jl")
include("operator.jl")
+include("image.jl")
function clip_line(r0, r1, l0, l1)
d_l = -(l1[1] - l0[1])
diff --git a/src/function.jl b/src/function.jl
index dfc7cfcd5f8cb98140cabe711f118d90b83d6272..6216bd9dc551dd62dcf12b23c3ad36e30f42ac07 100644
--- a/src/function.jl
+++ b/src/function.jl
@@ -257,31 +257,13 @@ function evaluate(df::Derivative, x)
end
-# TODO: inherit from some abstract function type
-struct ImageFunction{Img}
- mesh::Mesh
- img::Img
- cell::Base.RefValue{Int}
-end
-
-ImageFunction(mesh, img) =
- ImageFunction(mesh, img, Ref(1))
-
-bind!(f::ImageFunction, cell) = f.cell[] = cell
-img_coord(img, x) = (size(img, 1) - x[2] + 1, x[1])
-# TODO: precompute the offset from minimum mesh vertex
-evaluate(f::ImageFunction, xloc) =
- interpolate_bilinear(f.img,
- img_coord(f.img, elmap(f.mesh, f.cell[])(xloc) .+ (0.5, 0.5)))
-
-
struct FacetDivergence{F}
f::F
cell::Base.RefValue{Int}
edgemap::Dict{NTuple{2, Int}, Vector{Int}}
end
-# TODO: incomplete!
+# TODO: unfinished!
function FacetDivergence(f::FeFunction)
f.space.element == DP0() || throw(ArgumentError("unimplemented"))
mesh = f.space.mesh
@@ -324,59 +306,6 @@ end
bind!(f::FacetDivergence, cell) = f.cell[] = cell
-# evaluate the function on a rectangular grid of coordinates starting at (0.5,
-# 0.5) with integer spacing, i.e. pixel center coordinates of a bounding image.
-# indexing as for images: first index is y downwards, second index is x
-# rightwards
-# TODO: get rid of this hack as soon as we use size=() for scalar functions
-function sample(f::FeFunction)
- out = _sample(f)
- return f.space.size == (1,) ?
- reshape(out, size(out)[2:end]) : out
-end
-
-function _sample(f::FeFunction)
- # assumes dual grid
- mesh = f.space.mesh
-
- m0 = NTuple{2}(minimum(mesh.vertices, dims = 2))
- m1 = NTuple{2}(maximum(mesh.vertices, dims = 2))
-
- fcolon = ntuple(_ -> :, length(f.space.size))
-
- outsize = (f.space.size..., round.(Int, m1 .- m0 .+ 1)...)
- out = similar(f.data, outsize)
- for cell in cells(mesh)
- bind!(f, cell)
- A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
- view(mesh.vertices, :, view(mesh.cells, :, cell)))
-
- # TODO: cleanup when minimum/maximum with dims on StaticArrays becomes
- # type-stable
- I0 = round.(Int, SVector(minimum(A[1, :]), minimum(A[2, :])) .- m0) .+ 1
- I1 = round.(Int, SVector(maximum(A[1, :]), maximum(A[2, :])) .- m0) .+ 1
-
- for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2])
- x = SVector(Tuple(I) .- 1 .+ m0)
- # TODO: move to global-to-local function or inside-triangle test
- xloc = (A[:, SUnitRange(2, 3)] .- A[:, 1])::SArray \
- (x .- A[:, 1])
- ε = eps()
- if all(xloc .>= -ε) && sum(xloc) <= 1 + ε
- out[fcolon..., I] .= evaluate(f, xloc)
- end
- end
- end
-
- # convert to image indexing
- d = length(f.space.size)
- reverse!(out, dims = d + 2)
- out2 = permutedims(out, (ntuple(identity, d)..., d + 2, d + 1))
-
- return out2
-end
-
-
prolong!(new_f, old_cell, new_cells) =
prolong!(new_f, old_cell, new_cells, new_f.space.element)
diff --git a/src/image.jl b/src/image.jl
index 0ea64e5d2601d39f2826e0a3a89449d8a88e85fa..dfab9b856d2072fa3deef2b0cb460ee6d5ab4ecd 100644
--- a/src/image.jl
+++ b/src/image.jl
@@ -1,5 +1,27 @@
export interpolate_bilinear, halve, warp_backwards
+# ImageFunction
+
+# TODO: inherit from some abstract function type
+struct ImageFunction{Img}
+ mesh::Mesh
+ img::Img
+ cell::Base.RefValue{Int}
+end
+
+ImageFunction(mesh, img) =
+ ImageFunction(mesh, img, Ref(1))
+
+bind!(f::ImageFunction, cell) = f.cell[] = cell
+# transform coordinates to image/matrix indexing space
+img_coord(img, x) = (size(img, 1) - x[2] + 1, x[1])
+evaluate(f::ImageFunction, xloc) =
+ interpolate_bilinear(f.img,
+ img_coord(f.img, elmap(f.mesh, f.cell[])(xloc)))
+
+# 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
@@ -44,3 +66,311 @@ function halve(img)
return res
end
+
+# Sampling
+
+# evaluate the function on a rectangular grid of coordinates with integer
+# spacing bounded by rounded minimum and maximum vertex coordinates.
+# i.e. pixel center coordinates of an image. output indexing as for images:
+# first index is y downwards, second index is x rightwards
+# TODO: get rid of this hack as soon as we use size=() for scalar functions
+function sample(f::FeFunction)
+ out = _sample(f)
+ return f.space.size == (1,) ?
+ reshape(out, size(out)[2:end]) : out
+end
+
+function _sample(f::FeFunction)
+ # assumes dual grid
+ mesh = f.space.mesh
+
+ m0 = NTuple{2}(minimum(mesh.vertices, dims = 2))
+ m1 = NTuple{2}(maximum(mesh.vertices, dims = 2))
+
+ fcolon = ntuple(_ -> :, length(f.space.size))
+
+ outsize = (f.space.size..., round.(Int, m1 .- m0 .+ 1)...)
+ out = similar(f.data, outsize)
+ for cell in cells(mesh)
+ bind!(f, cell)
+ A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
+ view(mesh.vertices, :, view(mesh.cells, :, cell)))
+
+ # TODO: cleanup when minimum/maximum with dims on StaticArrays becomes
+ # type-stable
+ I0 = round.(Int, SVector(minimum(A[1, :]), minimum(A[2, :])) .- m0) .+ 1
+ I1 = round.(Int, SVector(maximum(A[1, :]), maximum(A[2, :])) .- m0) .+ 1
+
+ for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2])
+ x = SVector(Tuple(I) .- 1 .+ m0)
+ # TODO: move to global-to-local function or inside-triangle test
+ xloc = (A[:, SUnitRange(2, 3)] .- A[:, 1])::SArray \
+ (x .- A[:, 1])
+ ε = eps()
+ if all(xloc .>= -ε) && sum(xloc) <= 1 + ε
+ out[fcolon..., I] .= evaluate(f, xloc)
+ end
+ end
+ end
+
+ # convert to image indexing
+ d = length(f.space.size)
+ reverse!(out, dims = d + 2) # reverse y
+ out2 = permutedims(out, (ntuple(identity, d)..., d + 2, d + 1)) # flip xy
+
+ return out2
+end
+
+# Interpolation
+
+"""
+interprets `img` as a bilinearly interpolated continuous function and applies
+the default interpolation function for the discrete function u.
+"""
+function interpolate!(u, img)
+ f = ImageFunction(u.space.mesh, img)
+ interpolate!(u, @inline (x; f) -> f; f)
+end
+
+project_img(space::FeSpace, img) =
+ (u = FeFunction(space); project_img!(u, img))
+
+function project_img!(u::FeFunction, img)
+ d = 2 # domain dimension
+ space = u.space
+ mesh = space.mesh
+ f = ImageFunction(mesh, img)
+ opparams = (; f)
+
+ nrdims = prod(space.size)
+ # number of element dofs (i.e. local dofs not counting range dimensions)
+ nldofs = ndofs(space.element)
+
+ a(xloc, u, du, v, dv; f) = dot(u, v)
+ l(xloc, v, dv; f) = dot(f, v)
+
+ I = Float64[]
+ J = Float64[]
+ V = Float64[]
+ b = zeros(ndofs(space))
+
+ # mesh cells
+ for cell in cells(mesh)
+ foreach(f -> bind!(f, cell), opparams)
+ # cell map is assumed to be constant per cell
+ delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+ delmapinv = inv(delmap)
+ intel = abs(det(delmap))
+
+ p = ceil(Int, diam(mesh, cell))
+ qw, qx = quadrature_composite_lagrange_midpoint(p)
+ nqpts = length(qw) # number of quadrature points
+
+ qphi = zeros(nrdims, nrdims, nldofs, nqpts)
+ dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts)
+ for k in axes(qx, 2)
+ for r in 1:nrdims
+ qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k)))
+ dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k))))
+ end
+ end
+
+ # quadrature points
+ for k in axes(qx, 2)
+ xhat = SVector{d}(view(qx, :, k))
+ x = elmap(mesh, cell)(xhat)
+ opvalues = map(f -> evaluate(f, xhat), opparams)
+
+ # local test-function dofs
+ for jdim in 1:nrdims, ldofj in 1:nldofs
+ gdofj = space.dofmap[jdim, ldofj, cell]
+
+ phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj, k))
+ dphij = SArray{Tuple{space.size..., d}}(
+ SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj, k)) * delmapinv)
+
+ lv = qw[k] * l(x, phij, dphij; opvalues...) * intel
+ b[gdofj] += lv
+
+ # local trial-function dofs
+ for idim in 1:nrdims, ldofi in 1:nldofs
+ gdofi = space.dofmap[idim, ldofi, cell]
+
+ phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi, k))
+ dphii = SArray{Tuple{space.size..., d}}(
+ SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi, k)) * delmapinv)
+
+ av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel
+ push!(I, gdofi)
+ push!(J, gdofj)
+ push!(V, av)
+ end
+ end
+ end
+ end
+
+ ngdofs = ndofs(space)
+ A = sparse(I, J, V, ngdofs, ngdofs)
+
+ u.data .= A \ b
+ return u
+end
+
+project_img2(space::FeSpace, img) =
+ (u = FeFunction(space); project_img2!(u, img))
+
+# L1-stable quasi-interpolation operator
+# (https://arxiv.org/pdf/1505.06931.pdf)
+#
+# FIXME: currently only approximate quadrature by sampling on lagrange lattice
+# based on element size. Exact evaluation is tricky to implement (lots of
+# intersection handling).
+function project_img2!(u::FeFunction, img)
+ u.space.element == P1() ||
+ throw(ArgumentError("element unsupported"))
+ d = 2 # domain dimension
+ space = u.space
+ mesh = space.mesh
+ f = ImageFunction(mesh, img)
+ # count contributions to respective dof for subsequent averaging
+ gdofcount = zeros(Int, size(u.data))
+ area_refel = 0.5
+
+ u.data .= zero(eltype(u.data))
+ for cell in cells(mesh)
+ bind!(f, cell)
+
+ # size: nrdims × nldofs
+ # TODO: should be statically sized
+ gdofs = u.space.dofmap[:, :, cell]
+
+ p = ceil(Int, diam(mesh, cell))
+ qw, qx = quadrature_composite_lagrange_midpoint(p)
+ for k in axes(qx, 2)
+ xhat = SVector{d}(view(qx, :, k))
+ x = elmap(mesh, cell)(xhat)
+
+ # size: nrdims
+ value = evaluate(f, xhat)
+ # size: nldofs
+ phix = evaluate_basis(space.element, xhat)
+
+ # no integration element used
+ u.data[gdofs] .+= qw[k] ./ area_refel .* 3 .* value * phix'
+ #u.data[gdofs] .+= qw[k] ./ area_refel * value * (12 .* phix' .- 3)
+ end
+ gdofcount[gdofs] .+= 1
+ end
+ u.data ./= gdofcount
+ return u
+end
+
+# FIXME: unfinished / unused -> remove?
+function quadrature_cell_pixels(mesh, cell; m0)
+ d = ndims_domain(mesh)
+ A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
+ view(mesh.vertices, :, view(mesh.cells, :, cell)))
+
+ # TODO: cleanup when minimum/maximum with dims on StaticArrays becomes
+ # type-stable
+ I0 = round.(Int, SVector(minimum(A[1, :]), minimum(A[2, :])) .- m0) .+ 1
+ I1 = round.(Int, SVector(maximum(A[1, :]), maximum(A[2, :])) .- m0) .+ 1
+
+ #quadrature() = SA[1/6, 1/6, 1/6], SA[1/6 4/6 1/6; 1/6 1/6 4/6]
+ vertices = Float64[]
+ for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2])
+ x = SVector(Tuple(I) .- 1 .+ m0)
+ # TODO: move to global-to-local function or inside-triangle test
+ xloc = (A[:, SUnitRange(2, 3)] .- A[:, 1])::SArray \
+ (x .- A[:, 1])
+ ε = eps()
+ if all(xloc .>= -ε) && sum(xloc) <= 1 + ε
+ append!(vertices, xloc)
+ end
+ end
+ qw = ones(length(vertices) ÷ d) ./ length(vertices)
+ qx = reshape!(vertices, (d, :))
+end
+
+
+function interpolate_img_l2pixel!(u::FeFunction, img)
+ d = 2 # domain dimension
+ space = u.space
+ mesh = space.mesh
+ f = ImageFunction(mesh, img)
+ opparams = (; f)
+
+ nrdims = prod(space.size)
+ # number of element dofs (i.e. local dofs not counting range dimensions)
+ nldofs = ndofs(space.element)
+
+ a(xloc, u, du, v, dv; f) = dot(u, v)
+ l(xloc, v, dv; f) = dot(f, v)
+
+ I = Float64[]
+ J = Float64[]
+ V = Float64[]
+ b = zeros(ndofs(space))
+
+ # mesh cells
+ for cell in cells(mesh)
+ foreach(f -> bind!(f, cell), opparams)
+ # cell map is assumed to be constant per cell
+ delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+ delmapinv = inv(delmap)
+ intel = abs(det(delmap))
+
+ p = ceil(Int, diam(mesh, cell))
+ qw, qx = quadrature_composite_lagrange_midpoint(p)
+ nqpts = length(qw) # number of quadrature points
+
+ qphi = zeros(nrdims, nrdims, nldofs, nqpts)
+ dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts)
+ for k in axes(qx, 2)
+ for r in 1:nrdims
+ qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k)))
+ dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k))))
+ end
+ end
+
+ # quadrature points
+ for k in axes(qx, 2)
+ xhat = SVector{d}(view(qx, :, k))
+ x = elmap(mesh, cell)(xhat)
+ opvalues = map(f -> evaluate(f, xhat), opparams)
+
+ # local test-function dofs
+ for jdim in 1:nrdims, ldofj in 1:nldofs
+ gdofj = space.dofmap[jdim, ldofj, cell]
+
+ phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj, k))
+ dphij = SArray{Tuple{space.size..., d}}(
+ SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj, k)) * delmapinv)
+
+ lv = qw[k] * l(x, phij, dphij; opvalues...) * intel
+ b[gdofj] += lv
+
+ # local trial-function dofs
+ for idim in 1:nrdims, ldofi in 1:nldofs
+ gdofi = space.dofmap[idim, ldofi, cell]
+
+ phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi, k))
+ dphii = SArray{Tuple{space.size..., d}}(
+ SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi, k)) * delmapinv)
+
+ av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel
+ push!(I, gdofi)
+ push!(J, gdofj)
+ push!(V, av)
+ end
+ end
+ end
+ end
+
+ ngdofs = ndofs(space)
+ A = sparse(I, J, V, ngdofs, ngdofs)
+
+ u.data .= A \ b
+ return u
+
+end
diff --git a/src/mesh.jl b/src/mesh.jl
index a8c5c2e09056dee70516ff9d3836062e6cbc1e33..7177c01e4e8c679106e116dee01b946a365da014 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -1,5 +1,5 @@
export init_grid, init_hgrid, save, refine!, refine, cells, vertices,
- ndims_domain, ndims_space
+ ndims_domain, ndims_space, nvertices, nvertices_cell, ncells
using LinearAlgebra: norm
@@ -17,7 +17,9 @@ Base.show(io::IO, ::MIME"text/plain", x::HMesh) =
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)
diff --git a/src/operator.jl b/src/operator.jl
index c77afc464faedb674ec5a75ebffd70ca339dfe58..ea428e8ea6098d47ef335867dff539f83b139fa4 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -116,11 +116,6 @@ function assemble(space::FeSpace, a, l; params...)
return A, b
end
-from_img(arr) = permutedims(reverse(arr; dims = 1))
-
-project_img(space::FeSpace, img) =
- (u = FeFunction(space); project_img!(u, img))
-
# composite midpoint quadrature on 2d-lagrange point lattice
function quadrature_composite_lagrange_midpoint(p)
d_ = 2
@@ -139,133 +134,3 @@ function quadrature_composite_lagrange_midpoint(p)
return weights, points
end
-
-function project_img!(u::FeFunction, img)
- d = 2 # domain dimension
- space = u.space
- mesh = space.mesh
- f = ImageFunction(mesh, img)
- opparams = (; f)
-
- nrdims = prod(space.size)
- # number of element dofs (i.e. local dofs not counting range dimensions)
- nldofs = ndofs(space.element)
-
- a(xloc, u, du, v, dv; f) = dot(u, v)
- l(xloc, v, dv; f) = dot(f, v)
-
- I = Float64[]
- J = Float64[]
- V = Float64[]
- b = zeros(ndofs(space))
-
- # mesh cells
- for cell in cells(mesh)
- foreach(f -> bind!(f, cell), opparams)
- # cell map is assumed to be constant per cell
- delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
- delmapinv = inv(delmap)
- intel = abs(det(delmap))
-
- p = ceil(Int, diam(mesh, cell))
- qw, qx = quadrature_composite_lagrange_midpoint(p)
- nqpts = length(qw) # number of quadrature points
-
- qphi = zeros(nrdims, nrdims, nldofs, nqpts)
- dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts)
- for k in axes(qx, 2)
- for r in 1:nrdims
- qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k)))
- dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k))))
- end
- end
-
- # quadrature points
- for k in axes(qx, 2)
- xhat = SVector{d}(view(qx, :, k))
- x = elmap(mesh, cell)(xhat)
- opvalues = map(f -> evaluate(f, xhat), opparams)
-
- # local test-function dofs
- for jdim in 1:nrdims, ldofj in 1:nldofs
- gdofj = space.dofmap[jdim, ldofj, cell]
-
- phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj, k))
- dphij = SArray{Tuple{space.size..., d}}(
- SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj, k)) * delmapinv)
-
- lv = qw[k] * l(x, phij, dphij; opvalues...) * intel
- b[gdofj] += lv
-
- # local trial-function dofs
- for idim in 1:nrdims, ldofi in 1:nldofs
- gdofi = space.dofmap[idim, ldofi, cell]
-
- phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi, k))
- dphii = SArray{Tuple{space.size..., d}}(
- SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi, k)) * delmapinv)
-
- av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel
- push!(I, gdofi)
- push!(J, gdofj)
- push!(V, av)
- end
- end
- end
- end
-
- ngdofs = ndofs(space)
- A = sparse(I, J, V, ngdofs, ngdofs)
-
- u.data .= A \ b
- return u
-end
-
-project_img2(space::FeSpace, img) =
- (u = FeFunction(space); project_img2!(u, img))
-
-# L1-stable quasi-interpolation operator
-# (https://arxiv.org/pdf/1505.06931.pdf)
-#
-# FIXME: currently only approximate quadrature by sampling on lagrange lattice
-# based on element size. Exact evaluation is tricky to implement (lots of
-# intersection handling).
-function project_img2!(u::FeFunction, img)
- u.space.element == P1() ||
- throw(ArgumentError("element unsupported"))
- d = 2 # domain dimension
- space = u.space
- mesh = space.mesh
- f = ImageFunction(mesh, img)
- # count contributions to respective dof for subsequent averaging
- gdofcount = zeros(Int, size(u.data))
- area_refel = 0.5
-
- u.data .= zero(eltype(u.data))
- for cell in cells(mesh)
- bind!(f, cell)
-
- # size: nrdims × nldofs
- # TODO: should be statically sized
- gdofs = u.space.dofmap[:, :, cell]
-
- p = ceil(Int, diam(mesh, cell))
- qw, qx = quadrature_composite_lagrange_midpoint(p)
- for k in axes(qx, 2)
- xhat = SVector{d}(view(qx, :, k))
- x = elmap(mesh, cell)(xhat)
-
- # size: nrdims
- value = evaluate(f, xhat)
- # size: nldofs
- phix = evaluate_basis(space.element, xhat)
-
- # no integration element used
- u.data[gdofs] .+= qw[k] ./ area_refel .* 3 .* value * phix'
- #u.data[gdofs] .+= qw[k] ./ area_refel * value * (12 .* phix' .- 3)
- end
- gdofcount[gdofs] .+= 1
- end
- u.data ./= gdofcount
- return u
-end
diff --git a/test/Manifest.toml b/test/Manifest.toml
index 68a4dcd5eff907460e1c421cffd2020ad2e1aae5..25f1e8e3a1f6f3fd1e7233296bafddbf4dc6298f 100644
--- a/test/Manifest.toml
+++ b/test/Manifest.toml
@@ -1,26 +1,32 @@
# This file is machine-generated - editing it directly is not advised
-[[Base64]]
+julia_version = "1.7.0"
+manifest_format = "2.0"
+
+[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-[[InteractiveUtils]]
+[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-[[Logging]]
+[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-[[Markdown]]
+[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-[[Random]]
-deps = ["Serialization"]
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-[[Serialization]]
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-[[Test]]
+[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index 6bc28ce2831b0b5a34258575a4c49ca578c7e1f3..906ec15f8b29ea39d1416286715e93bda9e0de79 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,17 +12,19 @@ Random.seed!(0)
@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.]
- f_space = FeSpace(mesh, P1(), (1,))
- f = FeFunction(f_space)
+ u_space = FeSpace(mesh, P1(), (1,))
+ u = FeFunction(u_space)
- interpolate!(f, x -> interpolate_bilinear(img, x))
- img_sampled = sample(f)
+ interpolate!(u, img)
+ img_sampled = sample(u)
@test img == img_sampled
end