Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Interpolation between spaces exploiting tensor product structure #356

Open
IgorBaratta opened this issue Dec 3, 2021 · 1 comment
Open
Labels
enhancement New feature or request performance Performance related issues

Comments

@IgorBaratta
Copy link
Member

IgorBaratta commented Dec 3, 2021

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.

@IgorBaratta IgorBaratta added enhancement New feature or request performance Performance related issues labels Dec 3, 2021
@wence-
Copy link
Contributor

wence- commented Jan 13, 2022

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request performance Performance related issues
Projects
None yet
Development

No branches or pull requests

2 participants