code for obtaining subcell fluxes with desired sparsity pattern

[HennesHajduk]

As discussed, I've added code for obtaining subcell fluxes from an element vector that sums to zero. Note that there is an option bool UseDiagonalNbrs. If this is set to false, the dense stencil of a Q_1 element reduces to a cross stencil, similar to collocated methods. If true, the stencil of Q_1 remains unchanged, in either case the stencil is generally sparse.

Since I wanted to have this option and the integrals to be computed are non-standard, I found it easier to hard-code the entries, after tedious pen-and-paper calculations.

I'm performing the "assembly" on the reference element, and omitting the multiplication and subsequent division with the transformation constant. I strongly believe that this can also be done for severely curved elements, because Manuel observed that such integration errors are way less important than obtaining the right sparsity pattern.

The new functions are called from the main routine to illustrate how to use them. I haven't worked with Remhos in a while, so somebody else needs to apply this code where it's needed.