Algebraic simplification for the multipoint flux/stress approximation

[OmarDuran]

To eliminate the fluxes in the multipoint approximation, a block diagonal operator grad is constructed. The current implementation divides row and column permutation matrices into two linear applications, rows2blk_diag and cols2blk_diag, such that grad = rows2blk_diag @ grad_star @ cols2blk_diag, where grad_star is an operator being permuted. Performing the multipoint flux approximation with a direct construction of operator grad and possibly eliminating other global arrays will save computational memory and CPU time. This refactoring implies a deep modification of the files mpsa.py and mpfa.py. This issue reports that such an algebraic simplification may be necessary, particularly for larger problems.