add method for computing local approximation error

src/image.jl

@@ -443,6 +443,66 @@ function project_l2_pixel!(u::FeFunction, img)
     return u
 end
 
+"""
+    compute_local_error(u::FeFunction, img)
+
+Compute for every cell the squared L2-error between the mesh function `u` and
+the image `img`.
+Integrals are evaluated with a Lagrange quadrature rule with order dependent on
+the current cell diameter. Image point evaluations use bilinear interpolations.
+"""
+function compute_local_error(u::FeFunction, img)
+    d = 2
+    space = u.space
+    mesh = space.mesh
+
+    nrdims = prod(space.size)
+    # number of element dofs (i.e. local dofs not counting codomain dimensions)
+    nldofs = ndofs(space.element)
+
+    space_indicator = FeSpace(mesh, DP0(), (1,))
+    err = FeFunction(space_indicator; name = "err")
+
+    for cell in cells(mesh)
+        bind!(u, cell)
+        bind!(err, cell)
+        # cell map is assumed to be constant per cell
+        delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+        delmapinv = inv(delmap)
+        intel = abs(det(delmap))
+
+        p = ceil(Int, diam(mesh, cell))
+        qw, qx = quadrature_composite_lagrange_midpoint(p)
+        nqpts = length(qw) # number of quadrature points
+
+        qphi = zeros(nrdims, nrdims, nldofs, nqpts)
+        dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts)
+        for k in axes(qx, 2)
+            for r in 1:nrdims
+                qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k)))
+                dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k))))
+            end
+        end
+
+        # quadrature points
+        for k in axes(qx, 2)
+            xhat = SVector{d}(view(qx, :, k))
+            x = elmap(mesh, cell)(xhat)
+
+            # TODO: support more dimensions, evaluate_bilinear needs rework
+            u_val = first(evaluate(u, xhat))
+            img_val = evaluate_bilinear(img, x)
+
+            err_val = qw[k] * (u_val - img_val)^2 * intel
+
+            gdofs = space.dofmap[:, 1, cell]
+            err.data[gdofs] .= err_val
+        end
+
+    end
+    return err
+end
+
 function save_csv(path, img::AbstractArray{<:Any,2})
     df = DataFrame()
     for ci in CartesianIndices(img)