Infrastructure for matrix-valued problem

PR #392 should have the batteries for DG - but just for scalar problems. If we have vector valued problems, then our flux space becomes matrix-valued (think about displacement and stress). For convenience we could provide some handlers for handling matrix valued problems. This way DG can be implemented in a cleaner way than reformulating everything in Voigt notation beforehand.

Example reference: https://link.springer.com/content/pdf/10.1007/s10665-022-10238-3.pdf