rework image projection

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)
 
+        # size: nrdims × nldofs
+        # TODO: should be statically sized
+        gdofs = u.space.dofmap[:, :, 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)
@@ -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