From 7976b7f2907be4c68b7d2332083b9b718f2a8be6 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Mon, 28 Feb 2022 22:31:16 +0100
Subject: [PATCH] fixup residual error estimator
---
scripts/run_experiments.jl | 33 ++++++++++++++++++++-------------
1 file changed, 20 insertions(+), 13 deletions(-)
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 6e83010..a1c117d 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -616,8 +616,7 @@ function estimate_res!(st::L1L2TVState)
norm2(v) = dot(v, v)
cellf(hcell; g, u, nablau, p1, p2, tdata) =
- st.alpha2 *
- st.T(adjoint, tdata, st.alpha2 * st.T(tdata, u) - g) +
+ st.alpha2 * st.T(adjoint, tdata, st.T(tdata, u) - g) +
st.T(adjoint, tdata, p1) +
(isSnabla ? zero(u) : st.beta * u)
facetf(hfacet, n; g, u, nablau, p1, p2, tdata) =
@@ -645,21 +644,28 @@ function estimate_res!(st::L1L2TVState)
# cell contribution
hcell = diam(mesh, cell)
- delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+ cellmap = elmap(mesh, cell)
+ delmap = jacobian(cellmap, SA[0., 0.])
intel = abs(det(delmap))
- centroid = SArray{Tuple{ndims_domain(mesh)}}(mean(A, dims = 2))
- lcentroid = SA[1/3, 1/3]
-
- opvalues = map(f -> evaluate(f, lcentroid), fs)
-
- cellres = (isSnabla ? hcell^2 : hcell) *
- norm2(cellf(hcell; opvalues...)) .* intel
+ qw, qx = SemiSmoothNewton.quadrature()
+ val = 0.
+ for k in axes(qx, 2)
+ xhat = qx[:, k]
+ x = cellmap(xhat)
+ opvalues = map(f -> evaluate(f, xhat), fs)
+ val += qw[k] * intel *
+ (isSnabla ? hcell^2 : hcell) * norm2(cellf(hcell; opvalues...))
+ end
gdofs = space.dofmap[:, 1, cell]
- st.est.data[gdofs] .= cellres
+ st.est.data[gdofs] .= val
# facet contributions
+ centroid = SArray{Tuple{ndims_domain(mesh)}}(mean(A, dims = 2))
+ lcentroid = SA[1/3, 1/3]
+
+ opvalues = map(f -> evaluate(f, lcentroid), fs)
cross(v) = SA[-v[2], v[1]]
for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
fi = mesh.cells[i, cell]
@@ -693,7 +699,7 @@ function estimate_res!(st::L1L2TVState)
fj = mesh.cells[j, cell]
gdofs = space.dofmap[:, 1, cell]
- st.est.data[gdofs] .+= (isSnabla ? hfacet : inv(hfacet)) *
+ st.est.data[gdofs] .+= (isSnabla ? hfacet : inv(hfacet)) *
norm2(facetv[(fi, fj)]) * intel
end
end
@@ -849,13 +855,14 @@ function denoise(ctx)
# plot
i += 1
display(plot(grayclamp.(to_img(sample(st.u)))))
+ estimate_res!(st)
pvd[i] = save_step(i)
# refinement stop criterion
k_refine += 1
k_refine > n_refine && break
println("refine ...")
- estimate_res!(st)
+ #estimate_res!(st)
marked_cells = mark(st; theta = 0.5)
#marked_cells = Set(axes(mesh.cells, 2))
mesh, fs = refine(mesh, marked_cells;
--
GitLab