fix estimator per cell scaling

src/function.jl

@@ -156,6 +156,35 @@ end
 interpolate!(dst::FeFunction, ::DP1, expr::Function; params...) =
     interpolate!(dst, P1(), expr; params...)
 
+
+project!(dst::FeFunction, expr::Function; params...) =
+    project!(dst, dst.space.element, expr; params...)
+
+function project!(dst::FeFunction, ::DP0, expr; params...)
+    # same as interpolate!, but scaling by cell size
+    params = NamedTuple(params)
+    space = dst.space
+    mesh = space.mesh
+    for cell in cells(mesh)
+	delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+	intel = abs(det(delmap))
+
+	for f in params
+	    bind!(f, cell)
+	end
+
+	vertices = mesh.vertices[:, mesh.cells[:, cell]]
+	centroid = SArray{Tuple{ndims_domain(mesh)}}(mean(vertices, dims = 2))
+	lcentroid = SA[1/3, 1/3]
+
+	opvalues = map(f -> evaluate(f, lcentroid), params)
+
+	gdofs = space.dofmap[:, 1, cell]
+	dst.data[gdofs] .= myvec(expr(centroid; opvalues...)) .* intel
+    end
+
+end
+
 function bind!(f::FeFunction, cell)
     f.ldata .= vec(f.data[f.space.dofmap[:, :, cell]])
     return f

src/run.jl

@@ -200,7 +200,7 @@ function estimate!(ctx::L1L2TVContext)
     w = FeFunction(ctx.u.space)
     nablaw = nabla(w)
     solve_primal!(w, ctx)
-    interpolate!(ctx.est, estf; ctx.g, ctx.u, ctx.p1, ctx.p2, ctx.nablau, w, nablaw, ctx.tdata)
+    project!(ctx.est, estf; ctx.g, ctx.u, ctx.p1, ctx.p2, ctx.nablau, w, nablaw, ctx.tdata)
 end
 
 function refine(ctx::L1L2TVContext, marked_cells)