fixup residual error estimator

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]
@@ -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;