From 84f01d3172536d8bd2f1aef2fb2f546d91f0f767 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Wed, 14 Jul 2021 12:44:35 +0200
Subject: [PATCH] fix estimator per cell scaling
---
src/function.jl | 29 +++++++++++++++++++++++++++++
src/run.jl | 2 +-
2 files changed, 30 insertions(+), 1 deletion(-)
diff --git a/src/function.jl b/src/function.jl
index 2ea97ed..f0221fe 100644
--- a/src/function.jl
+++ b/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
diff --git a/src/run.jl b/src/run.jl
index cc2bbfd..ef2083d 100644
--- a/src/run.jl
+++ b/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)
--
GitLab