Skip to content
Snippets Groups Projects
Commit 20736d4b authored by Stephan Hilb's avatar Stephan Hilb
Browse files

implement l1-stable quasi-interpolation

parent f31f6cb3
No related branches found
No related tags found
No related merge requests found
......@@ -3,7 +3,7 @@ using LinearAlgebra: det, dot
using StaticArrays: SA, SArray, MArray
using ForwardDiff: jacobian
export Poisson, L2Projection, init_point!, assemble, project_img
export Poisson, L2Projection, init_point!, assemble, project_img, project_img2
abstract type Operator end
......@@ -119,22 +119,8 @@ 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)
# composite midpoint quadrature on lagrange point lattice
function quadrature(p)
function quadrature_composite_lagrange_midpoint(p)
d_ = 2
n = binomial(p + 2, 2)
weights = Vector{Float64}(undef, n)
......@@ -152,11 +138,24 @@ function project_img!(u::FeFunction, img)
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))
gdof = LinearIndices((nrdims, ndofs(space)))
# mesh cells
for cell in cells(mesh)
......@@ -167,7 +166,7 @@ function project_img!(u::FeFunction, img)
intel = abs(det(delmap))
p = ceil(Int, diam(mesh, cell))
qw, qx = quadrature(p)
qw, qx = quadrature_composite_lagrange_midpoint(p)
nqpts = length(qw) # number of quadrature points
qphi = zeros(nrdims, nrdims, nldofs, nqpts)
......@@ -219,3 +218,56 @@ function project_img!(u::FeFunction, img)
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)
nrdims = prod(space.size)
# number of element dofs (i.e. local dofs not counting range dimensions)
nldofs = ndofs(space.element)
# count contributions to respective dof for subsequent averaging
gdofcount = zeros(Int, size(u.data))
# mesh cells
for cell in cells(mesh)
bind!(f, cell)
p = ceil(Int, diam(mesh, cell))
# we actually only use the lagrange lattice points
qw, qx = quadrature_composite_lagrange_midpoint(p)
nqpts = length(qw) # number of quadrature points
qphi = zeros(nldofs, nqpts)
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, SVector{d}(view(qx, :, k)))
# size: nrdims × nldofs
gdofs = u.space.dofmap[:, :, cell]
u.data[gdofs] .+= value * (12 .* phix' .- 3)
gdofcount[gdofs] .+= 1
end
end
u.data ./= gdofcount
return u
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment