residual calculation needs to be done when boundary degrees of freedom have not been removed (since their removal takes place in the Mult step in order to introduce the proper restoring force). Whereas jacobian calculation needs to always work from the boundary aware operators (these are sparse matrices)
The other thing to note here is that this seems contrary to a matrix free design. While true it is common practive to assemble linear terms since they are cheap. We still never assemble the non linear matrix form.
In order to maintain memory efficienty I have implimented a matrix free SchurComplement operator as well as an operator which uses a few iterations of GMRES to approxinate the inverse of any general operator.
A custom integrator is required to handle the theta^n term in the lane emden equation, that is written as NonlinearPowerIntegrator which is a mfem::NonlinearFormIntegrator and defines methods to assemble its element vector (function value) and element gradient matrix (jacobian). This is then, along with built in mfem vectors for M Q and D, incorporated into the PolytropeOperator which defines methods for Mult (calculate the residuals of the variational form) and GetGradient (find the jacobian of the system)
