Fix fieldsplit with Cofunction right hand side.

[JHopeCollins]

Problem description

As described in #3368, it currently isn't possible to use PCFieldsplit if the form contains a Cofunction.
This is because the way we split the form is by traversing the dag with a ufl.MultiFunction and splitting each node as it is visited. However, none of the current handlers will deal with a Cofunction.

This PR adds a Cofunction handler to the ExtractSubBlock multifunction to do this splitting.
Fixes #3368

What does the Cofunction handler need to achieve?

A fieldsplit solver needs the right hand side from the corresponding field. For a Cofunction this means the values in the corresponding subfunction. So when we split a form the return value needs to live in the correct subspace, and point to the numerical values from the full space.

How is this implemented?

If all terms in the form have 1 argument (i.e. a FormSum of Cofunctions and 1-forms), then life is simple. We just check which field(s) is being requested, and make a new Cofunction in that (possibly mixed) subspace using the Dats from the original Cofunction.

If some terms in the full Form have 2 arguments, then life is a little more complicated, and we have two possibilities.

1. If we're being asked for a block on the diagonal, then a) numerically this is probably going to be inverted and so needs a right hand side and b) mathematically this acts on a member of a function space and returns an assembled thing in the dual space - this is exactly what the Cofunction is (I'm jetlagged and I've probably slaughtered the dual/primal space technicalities. Edits to fix the explanation welcome). In this case we do the same as for the 1-form case i.e. we just make a new smaller Cofunction that views the appropriate parts of the full Cofunction.

2. If we're being asked for an off-diagonal block, then a) numerically this doesn't make sense to invert and hence to have a right hand side, and b) mathematically the test and trial spaces are different, which doesn't make sense to have a corresponding Cofunction (I think?). In this case we just return a ZeroBaseForm which will get optimised out before we do any actual computational work.

Other changes

Because we can now split a Cofunction or FormSum, not just a Form, I had to widen the logic in split_form for checking if the form returned by ExtractSubBlock is empty or not. It works but it might not be the best/most robust way of checking.

What if this all breaks in the future?

There's a test that solves a mixed formulation of the wave equation using a Schur factorisation with both a 1-form right hand side and an equivalent Cofunction right hand side, and checks the results match close to machine precision.
The Schur complement is built by asking PETSc to do the brute force thing. Computationally suboptimal but the mesh is tiny, it reduces the amount of other "stuff" needed, and testing implementations is basically what the brute force thing is for.

[JHopeCollins]

@stephankramer I see that you were also having issues with this in thetisproject/thetis#376. If you/someone else from the Thetis team could checkout this branch and make sure it fixes things for you then that would be really helpful.

[github-actions]

	Tests	Passed ✅	Skipped ⏭️	Failed ❌
Firedrake real	8159 ran	7484 passed	675 skipped	0 failed

[github-actions]

	Tests	Passed ✅	Skipped ⏭️	Failed ❌
Firedrake complex	8153 ran	6679 passed	1474 skipped	0 failed

[connorjward]

LGTM

[JHopeCollins]

@pbrubeck unless you have any other issues I'd like to merge this PR.

all suggestions addressed