toarray assumes a knot multiplicity of 1

[jowezarek]

from    mpi4py              import MPI
import  numpy               as np
import  matplotlib.pyplot   as plt

from    sympde.expr                 import BilinearForm, integral
from    sympde.topology             import elements_of, Cube, Derham

from    psydac.api.discretization   import discretize
from    psydac.api.settings         import PSYDAC_BACKEND_GPYCCEL
from    psydac.linalg.basic         import LinearOperator
from    psydac.linalg.block         import BlockVectorSpace
from    psydac.linalg.stencil       import StencilVectorSpace

comm     = MPI.COMM_WORLD
backend  = PSYDAC_BACKEND_GPYCCEL

mult     = [2, 1, 1]

ncells   = [3, 2, 2]
degree   = [2, 1, 1]
periodic = [False]*3

domain   = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
derham   = Derham(domain)

domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm)
derham_h = discretize(derham, domain_h, degree=degree, multiplicity=mult)

u0, v0   = elements_of(derham.V0, names='u0, v0')

a0       = BilinearForm((u0, v0), integral(domain, u0*v0))
a0_h     = discretize(a0, domain_h, (derham_h.V0, derham_h.V0), backend=backend)
M0       = a0_h.assemble()

plt.spy(M0.toarray())
plt.show()
#plt.spy(toarray(M0))
#plt.show()

This code assembles the mass matrix corresponding to the first function space of a 3D de Rham sequence, where the knot sequence in x-diretcion is of multiplicity 2 instead of the default 1.
plt.spy(M0.toarray()) plots a non-symmetric matrix:

I wrote a custom (slow) toarray() function that recovers the (hopefully / probably) true symmetric structure:

The issue is most probably related to the fact that the pads value should at some point be multiplied with the knot multiplicity.