I have managed to sort of implement this using Gridap’s low-level API (as described in tutorial 13).
The implementation is available here: ee4375_fem_ta/High-Voltage Cable - Gridap.ipynb at main · gijswl/ee4375_fem_ta · GitHub
After defining the linear and bilinear terms of the weak form, we can use the low-level API to derive the cell contributions and assemble the linear system. This linear system can then be freely modified (e.g., extended with additional rows and columns). The cell contributions to these extra vectors can also be calculated with the low-level API in quite a nice way.