Add support for bilinear interpolation in the distributed remapper

[Sbozzolo]

The current non-conservative remapping we are using can change the extrama within an element. It would be good to support bilinear interpolation too, so that we can preserve the extrema.

The way the remapper currently works is this:

• For each target point, compute the two interpolation vector, one for each horizontal direction. These vectors have size Nq (the number of quadrature points). When you compute one interpolation vector per target point, you obtain an interpolation matrix.
• Also compute vertical interpolation matrices for linear interpolation
• Loop over all the points, compute something along the lines of V * W1 * W2 * F, with V vertical weights, W1/W2 horizontal weights, F field value.

To substitute bilinear interpolation, one would have to compute different interpolation weight matrices. These matrices have the challenge that not all the nodal points are relevant, so they wil require some seaching to build. It would be good to make them still with size Nq (but block matrices) so that both interpolation schemes can be used and the only change is defining the interpolation matrix.