From 57fb81e4b19d43c5b5f905228a6946d70c2b33a4 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Thu, 5 May 2022 14:27:12 +0200
Subject: [PATCH] add method for computing local approximation error
---
src/image.jl | 60 ++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 60 insertions(+)
diff --git a/src/image.jl b/src/image.jl
index eb44867..10e3740 100644
--- a/src/image.jl
+++ b/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)
--
GitLab