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

finish l2 pixel projection

parent c3e3716e
No related branches found
No related tags found
No related merge requests found
...@@ -12,7 +12,8 @@ using Plots ...@@ -12,7 +12,8 @@ using Plots
using SemiSmoothNewton using SemiSmoothNewton
using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs
using SemiSmoothNewton: project_l2_lagrange!, project_qi_lagrange!, project! using SemiSmoothNewton: project!, project_l2_lagrange!, project_qi_lagrange!,
project_l2_pixel!
using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save
include("util.jl") include("util.jl")
......
...@@ -295,6 +295,55 @@ function quadrature_cell_pixels(mesh, cell; m0) ...@@ -295,6 +295,55 @@ function quadrature_cell_pixels(mesh, cell; m0)
qx = reshape!(vertices, (d, :)) qx = reshape!(vertices, (d, :))
end end
struct PixelIterator{A, G}
"cell vertices"
A::A
"CartesianIndices covering the cell"
grid::G
end
Base.IteratorSize(::PixelIterator) = Base.SizeUnknown()
Base.eltype(it::PixelIterator) = eltype(it.grid)
"""
PixelIterator(mesh, cell)
PixelIterator(mesh, A)
Return an iterator over all integer cartesian indices (corresponding to pixel
centers) intersecting a given `cell` or geometry `A`.
"""
function PixelIterator(mesh::Mesh, cell::Int)
A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
view(mesh.vertices, :, view(mesh.cells, :, cell)))
return PixelIterator(mesh, A)
end
function PixelIterator(mesh::Mesh, A)
# rounding is probably correct since we don't miss any coordinates even
# when we round inwards to the element
I0 = round.(Int, SVector(minimum(A[1, :]), minimum(A[2, :])))
I1 = round.(Int, SVector(maximum(A[1, :]), maximum(A[2, :])))
return PixelIterator(A, CartesianIndex(I0...):CartesianIndex(I1...))
end
function _intersects(it::PixelIterator, I)
x = SVector(Tuple(I))
# TODO: move to global-to-local function or inside-triangle test
xloc = (it.A[:, SUnitRange(2, 3)] .- it.A[:, 1])::SArray \
(x .- it.A[:, 1])
ε = eps()
return all(xloc .>= -ε) && sum(xloc) <= 1 + ε
end
function Base.iterate(it::PixelIterator, state...)
y = iterate(it.grid, state...)
isnothing(y) && return nothing
while !_intersects(it, first(y))
y = iterate(it.grid, last(y))
isnothing(y) && return nothing
end
return y
end
function project_l2_pixel!(u::FeFunction, img) function project_l2_pixel!(u::FeFunction, img)
d = 2 # domain dimension d = 2 # domain dimension
...@@ -307,6 +356,20 @@ function project_l2_pixel!(u::FeFunction, img) ...@@ -307,6 +356,20 @@ function project_l2_pixel!(u::FeFunction, img)
# number of element dofs (i.e. local dofs not counting range dimensions) # number of element dofs (i.e. local dofs not counting range dimensions)
nldofs = ndofs(space.element) nldofs = ndofs(space.element)
# first loop to count number of triangles intersecting pixel evaluation
# points
ncells = zeros(axes(img))
for cell in cells(mesh)
pixels = PixelIterator(mesh, cell)
for I in pixels
ncells[I] += 1
end
end
all(!iszero, ncells) ||
throw(ErrorException("not all pixels are covered by a cell"))
pixelweights = inv.(ncells)
# typical L2 projection of f
a(xloc, u, du, v, dv; f) = dot(u, v) a(xloc, u, du, v, dv; f) = dot(u, v)
l(xloc, v, dv; f) = dot(f, v) l(xloc, v, dv; f) = dot(f, v)
...@@ -315,7 +378,10 @@ function project_l2_pixel!(u::FeFunction, img) ...@@ -315,7 +378,10 @@ function project_l2_pixel!(u::FeFunction, img)
V = Float64[] V = Float64[]
b = zeros(ndofs(space)) b = zeros(ndofs(space))
# mesh cells # buffer for basis value/derivative at a single point
qphi = zero(MArray{Tuple{nrdims, nrdims, nldofs}})
dqphi = zero(MArray{Tuple{nrdims, d, nrdims, nldofs}})
for cell in cells(mesh) for cell in cells(mesh)
foreach(f -> bind!(f, cell), opparams) foreach(f -> bind!(f, cell), opparams)
# cell map is assumed to be constant per cell # cell map is assumed to be constant per cell
...@@ -323,45 +389,49 @@ function project_l2_pixel!(u::FeFunction, img) ...@@ -323,45 +389,49 @@ function project_l2_pixel!(u::FeFunction, img)
delmapinv = inv(delmap) delmapinv = inv(delmap)
intel = abs(det(delmap)) intel = abs(det(delmap))
p = ceil(Int, diam(mesh, cell)) #p = ceil(Int, diam(mesh, cell))
qw, qx = quadrature_composite_lagrange_midpoint(p) #qw, qx = quadrature_composite_lagrange_midpoint(p)
nqpts = length(qw) # number of quadrature points #nqpts = length(qw) # number of quadrature points
qphi = zeros(nrdims, nrdims, nldofs, nqpts) # pixels as quadrature points
dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts) pixels = PixelIterator(mesh, cell)
for k in axes(qx, 2) for P in pixels
x = SVector(Tuple(P))
# TODO: move to global-to-local function or inside-triangle test
xhat = (pixels.A[:, SUnitRange(2, 3)] .- pixels.A[:, 1])::SArray \
(x .- pixels.A[:, 1])
opvalues = map(f -> evaluate(f, xhat), opparams)
# fill basis value/derivative buffer
for r in 1:nrdims for r in 1:nrdims
qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k))) qphi[r, r, :] .=
dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k)))) evaluate_basis(space.element, SVector{d}(xhat))
end dqphi[r, :, r, :]::MArray{Tuple{d, nldofs}, Float64, 2, d * nldofs} .=
transpose(jacobian(
x -> evaluate_basis(space.element, x),
SVector{d}(xhat)))
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 # local test-function dofs
for jdim in 1:nrdims, ldofj in 1:nldofs for jdim in 1:nrdims, ldofj in 1:nldofs
gdofj = space.dofmap[jdim, ldofj, cell] gdofj = space.dofmap[jdim, ldofj, cell]
phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj, k)) phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj))
dphij = SArray{Tuple{space.size..., d}}( dphij = SArray{Tuple{space.size..., d}}(
SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj, k)) * delmapinv) SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj)) * delmapinv)
lv = qw[k] * l(x, phij, dphij; opvalues...) * intel lv = pixelweights[P] * l(x, phij, dphij; opvalues...) * intel
b[gdofj] += lv b[gdofj] += lv
# local trial-function dofs # local trial-function dofs
for idim in 1:nrdims, ldofi in 1:nldofs for idim in 1:nrdims, ldofi in 1:nldofs
gdofi = space.dofmap[idim, ldofi, cell] gdofi = space.dofmap[idim, ldofi, cell]
phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi, k)) phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi))
dphii = SArray{Tuple{space.size..., d}}( dphii = SArray{Tuple{space.size..., d}}(
SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi, k)) * delmapinv) SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi)) * delmapinv)
av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel av = pixelweights[P] * a(x, phii, dphii, phij, dphij; opvalues...) * intel
push!(I, gdofi) push!(I, gdofi)
push!(J, gdofj) push!(J, gdofj)
push!(V, av) push!(V, av)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment