[firedrake] Assembling Jacobian within TAO

Hi all, So I am attempting to use TAO's mixed complementarity solvers and it basically requires a formGradient/Constraints routine and a formJacobian routine. I am solving a variational inequality via TAO and want the linear operator A to be my jacobian. However, I don't think it was correctly done because my semi smooth linear solver does not converge (but it should!). Attached is the code. Line 108 is what TaoSetJacobianRoutine() calls. Basically I want to assemble the linear operator into the matrix petsc_J. Is what I have the correct implementation? run the code as:

python MCP_poisson.py 100 100 2

Thanks, Justin