add unfinished p1_project! work

scripts/run_experiments.jl

@@ -151,13 +151,36 @@ function OptFlowState(mesh;
 	est, g, u, p1, p2, du, dp1, dp2)
 end
 
-function p1_project!(p1, alpha1)
-    p1.space.element == DP0() || p1.space.element == P1() ||
-        p1.space.element == DP1() ||
-        throw(ArgumentError("element unsupported"))
+"calculate l1 norm on reference triangle from langrange dofs"
+function p1_ref_l1_expr(ldofs)
+    u0, u1, u2 = ldofs
+    return (
+            (u0 - u1) * abs(u2) * u2 ^ 2 -
+            ((u0 - u2) * abs(u1) * u1 ^ 2 -
+                (u1 - u2) * abs(u0) * u0 ^ 2)) /
+        (6 * ((u0 ^ 2 + u1 * u2) * (u1 - u2) - (u1 ^ 2 - u2 ^ 2) * u0))
+end
+
+p1_project!(p1, alpha1) = p1_project!(p1, alpha1, p1.space.element)
+p1_project!(p1, alpha1, _) = throw(ArgumentError("element unsupported"))
+
+# FIXME: probably only correct for DP0?
+# TODO: finish correct discrete projection?!
+function p1_project!(p1, alpha1, ::Union{DP0, P1, DP1})
     p1.data .= clamp.(p1.data, -alpha1, alpha1)
 end
 
+#function p1_project!(p1, alpha1, ::Union{P1, DP1})
+#    mesh = p1.space.mesh
+#    for cell in cells(mesh)
+#        bind!(p1, cell)
+#        F = elmap(mesh, cell)
+#
+#        ForwardDiff.gradient(p1_ref_l1_expr, p1.ldata)
+#
+#    end
+#end
+
 function p2_project!(p2, lambda)
     p2.space.element == DP0() || p2.space.element == DP1() ||
         p2.space.element == P1() ||
@@ -251,6 +274,7 @@ function step!(ctx::L1L2TVState)
     ctx.p2.data .+= theta * ctx.dp2.data
 
     # reproject p1, p2
+    # FIXME: the p1 projection destroys the primal reconstruction
     p1_project!(ctx.p1, ctx.alpha1)
     p2_project!(ctx.p2, ctx.lambda)
 end