diff --git a/src/operator.jl b/src/operator.jl index 87c3710161f0f14c0e2786b8932f30a35eb500b9..c77afc464faedb674ec5a75ebffd70ca339dfe58 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