Assembling matrix via inner(div(u), div(v)) failed about shenfun HOT 4 CLOSED

from sympy import symbols, cos, sin
from mpi4py import MPI
from shenfun import inner, grad, div, TestFunction, TrialFunction, Array,\
    Basis, TensorProductSpace, VectorTensorProductSpace, BlockMatrix
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD

family = 'legendre'

# Size of discretization
n = 40
N = [n, n+1]

# Bases and spaces
BX_hom = Basis(N[0], family=family, bc=(0, 0))
BX_inhom = Basis(N[0], family=family, bc=(0, 0.1))
BY_hom = Basis(N[1], family=family, bc=(0, 0))

# Function spaces
subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T_hom = TensorProductSpace(subcomms, (BX_hom, BY_hom))
T_inhom = TensorProductSpace(subcomms, (BX_inhom, BY_hom))

V_hom = VectorTensorProductSpace([T_hom, T_hom])
V_inhom = VectorTensorProductSpace([T_inhom, T_inhom])
V_mixed = VectorTensorProductSpace([T_hom, T_inhom])

# Choose the BCs you want
# working
V = V_hom

# not working
# V = V_inhom
# V = V_mixed

# Functions
u = TrialFunction(V)
v = TestFunction(V)

# Compute matrices
# Assemble standard term
A = inner(grad(u), grad(v))

# Assemble (div(u), div(v))-term
# non-working
B = inner(div(u), div(v)) # raises an error

# Compute BlockMatrix
M = BlockMatrix(A + B)

# Compute right hand side
# Use sympy to compute a rhs
x, y = symbols("x,y")
fe = (sin(x**2)*(1-y)*(1+y), y**2*(cos(x-y**2)))
fj = Array(V, buffer=fe)

b = inner(v, fj)

# Compute solution
u_hat = M.solve(b) # This line fails because inconsistent shapes
u_sol = u_hat.backward()