From 68edb2a968b311bd1391473ab97cbb3b48b49ec8 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Sat, 27 Nov 2021 18:55:54 +0100
Subject: [PATCH] rework image projection
---
src/operator.jl | 29 ++++++++++++-----------------
1 file changed, 12 insertions(+), 17 deletions(-)
diff --git a/src/operator.jl b/src/operator.jl
index 87c3710..c77afc4 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -132,7 +132,7 @@ function quadrature_composite_lagrange_midpoint(p)
for I in Iterators.product(ntuple(_ -> 0:p, d_)...)
sum(Tuple(I)) > p && continue
k += 1
- weights[k] = 1 / n
+ weights[k] = 1 / n / factorial(d_)
points[1, k] = I[1] / p
points[2, k] = I[2] / p
end
@@ -237,25 +237,20 @@ function project_img2!(u::FeFunction, img)
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)
-
- u.data .= 0
# count contributions to respective dof for subsequent averaging
gdofcount = zeros(Int, size(u.data))
+ area_refel = 0.5
- # mesh cells
+ u.data .= zero(eltype(u.data))
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
+ # size: nrdims × nldofs
+ # TODO: should be statically sized
+ gdofs = u.space.dofmap[:, :, cell]
- qphi = zeros(nldofs, nqpts)
+ 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)
@@ -265,11 +260,11 @@ function project_img2!(u::FeFunction, img)
# size: nldofs
phix = evaluate_basis(space.element, xhat)
- # size: nrdims × nldofs
- gdofs = u.space.dofmap[:, :, cell]
- u.data[gdofs] .+= value * (12 .* phix' .- 3)
- gdofcount[gdofs] .+= 1
+ # 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
--
GitLab