io: fix saving/loading of HDiv/HCurl functions on a high-order mesh

[ksagiyam]

Fix #3781.

As @colinjcotter suggested on slack, we map HDiv/HCurl function defined on a high-order mesh to a DG function defined on an affine mesh using the contravariant/covariant Piola map.

Notes:

If we attach the Piola map to the vector DG element (using WithMapping), the DG DoF values will then depend on the reference-physical element map, which is inconvenient for checkpointing.

We need to have users pass a representation of the affine mesh explicitly as the mere CG1 function on a high-order mesh does not necessarily represent a sane affine mesh; e.g., a single CG2 mesh cell can look like the following, and CG1 mesh connecting vertices would collapse onto a line:

    .   --o --
 /                \
          o         .
|      /     \      |
x--o--x       x--o--x

We need to have users pass the quadrature degree used when projecting to the DG space.

[github-actions]

	Tests	Passed ✅	Skipped ⏭️	Failed ❌
Firedrake complex	8084 ran	6501 passed	1583 skipped	0 failed

[github-actions]

	Tests	Passed ✅	Skipped ⏭️	Failed ❌
Firedrake real	8090 ran	7304 passed	786 skipped	0 failed