Interpolation between spaces exploiting tensor product structure
Opened this issue · 1 comments
The number of floating point operations per cell is currently O((P0+1)^d *(P1+1)^d) =~ O(P^2*d)
where P0 is the degree of interpolation from and P1 the degree for interpolation to.
If the pair of elements have tensor product factorization, this number can be reduced to ~ P^d+1
, if only factor matrices are computed.
I'm not sure if basix should return this 1d matrices, or the user should request the interpolation matrix for a pair of factors instead.
FWIW, I think if you're going to hand back tabulations to support this you need to hand back a matrix that contains the structure (that is, it can't just be a dense array) that implements the matrix-vector product efficiently.
If you hand back 1d factors, it is now the caller's job to glue the matvec together in the right way, which I guess you would want to wrap up. But now that makes writing generic code a bit harder (because you have to dispatch on whether or not the element has tensor product structure).
I think the nice way to do this is either extend the interpolate interface to offer an FiniteElement::interpolate_into(const xtensor<double>& from, FiniteElement target, xtensor<double>& to)
that doesn't build the interpolation matrix but instead applies the the action of the interpolation matrix on the coefficients from
to produce coefficients to
in the target space.
You could type-specialise this to pairs of finite elements that support fast algorithms and fall back to the dense mat-vec for everything else (to handle everything there's a cartesian product explosion of cases which is awkward).
Alternately, the type of the return value of FiniteElement::interpolation_matrix
would be some abstract matrix type that implements mat-vec (and perhaps as_xtensor
or something). But now you have to implement the relevant algebra on these matrices for composition and so forth.
I think all of the above applies, mutatis mutandis, to FiniteElement::tabulate
as well.