Add function to fully fill a PETSc matrix with zeros

After preallocating the correct amount of space, we fill the entire
matrix with zeros so that on first assembly PETSc does not compact out
some of the space we said we'd need.
Fully fill matrices with zeros on construction

The sparsity pattern we specify is the right one, but PETSc compresses
the allocation on matrix assembly to the number of "used" nonzeros.
This can result in requiring reallocation if (for example) the values
change because we're actually assembling the linearisation of a
non-linear operator.  To avoid this, when we first build a matrix, fill
all the entries we might ever use with zeros.
Rework sparsity construction

Unify sparsity construction for the MPI and non-MPI cases.
Additionally, simplify the extruded sparsity construction code to reduce
repetition.  As an additional bonus, we now correctly build sparsity
patterns for all extruded cases in parallel as well as serial.

Performance is unchanged from the previous code, with some marginal
improvements in parallel extruded VFS cases (but it's basically a
wash).
Use BAIJ matrices where possible

When building a matrix where rdim == cdim > 1, we now use
block-sparse (BAIJ) matrices.  This reduces the size of the sparsity
pattern we compute by a factor of rdim*cdim and allows for more
efficient matrix storage and insertion.  As a consequence, we no longer
create sequential matrices with arrays (these are not supported for
PETSc BAIJ matrices).  The parallel and sequential code therefore
becomes more closely aligned.
Remove unused variable
Update tests in light of BAIJ changes

We no longer have any access to the matrix data array directly.
Error on device backends for unimplemented Mat behaviour

We drop support for matrices on OpenCL completely, and do not yet
implement block-sparse matrices on Cuda.  Tests modified appropriately
to skip tests on device backends where the functionality is
unsupported.