Kronecker products

[dlfivefifty]

I think kron should kronecker the indices as well. That is

x = Inclusion(0..1)
y = Inclusion(2..3)

xy = kron(x,y)
axes(xy) == (Inclusion((0..1) × (2..3)),) # a rectangle
eltype(xy) == 
xy[SVector(0.1,0.2)] == 0.1*0.2

This is inconsistent with kron for arrays, but I don't think there is any other sense in which this works. We may need to call it something else, say tensorproduct or ⊗ to allow access for standard arrays.

It also then matches the physical sense and we can do for example tensor product bases:

T = Chebyshev()
T2 = kron(T,T)
axes(T2) == (Inclusion((0..1)^2),(1:∞) ⊗ (1:∞))
T2[SVector(0.1,0.2), SVector(1,5)] == T[0.1,1] * T[0.2,5]