How to split an AbstractQuasiMatrix by coordinate, not basis function?

[jagot]

Related to #8

I started thinking again how to implement overlapping domains; when solving PDEs like the Schrödinger equation, different but equivalent formulations are more/less stiff in different parts of space. The formulations are related via a unitary transformation. One can then solve two problems in slightly overlapping domains and communicating the overlap between them every time step.

Setting up the equations is easy, as soon as the correct AbstractQuasiMatrices are known, the question is rather how to split the whole domain into subdomains. Imagine

R = FiniteDifferences(0.0..100, δx=1.0)

I now want to split R into subdomains R1 and R2, covering 0.0..30 and 20.0..100, respectively. I also want all materialized matrices to be of appropriate sizes, i.e. not full matrices with zeros outside their respective subdomains. This can be accomplished by

R1 = R[:,1:31]
R2 = R[:,21:end]

but then the first axes of R1 and R2 will be "wrong" and if R is not a uniform grid, as in this case, it is not so trivial to determine which basis functions need to be included. If on the other hand I do

R1 = R[Inclusion(0.0..30),:]
R2 = R[Inclusion(20.0..100),:]

then the first axes are correct but all the basis functions are still there, so the materialized matrices will still be too large.

• What would be a good interface for this?
• This obviously only makes sense for local bases, e.g. something like orthogonal polynomials would not fit into this. It's also not clear how this would translate to B-splines that have compact support, but are not orthogonal, i.e. the overlap matrix is banded.

[dlfivefifty]

I think you are trying to do the right thing but there's some bugs floating around where I assumed Inclusion was "everything", which is not the case.