algebraic operators do not handle sympy/Matrix

[ratnania]

The algebraic operators dot, cross, inner and outer do not handle correctly Matrix objects from sympy.

We consider the following

from sympde.topology import VectorFunctionSpace, Cube, element_of
from sympde.calculus import grad, dot, inner, outer, cross, div, Transpose

from sympy import Matrix

domain = Cube()

I = Matrix([[1, 2, 3], [0, 1, 0], [0, 0, 1]])

V = VectorFunctionSpace('V', domain)

u,v = [element_of(V, name=i) for i in ['u', 'v']]

Then all the following operations lead to the same error

>>> A = Matrix([[1, 2, 3], [2, 1, 0], [3, 0, 1]])

>>> dot(A*u, v)
...
AttributeError: 'MutableDenseMatrix' object has no attribute 'is_commutative'
>>> cross(A*u, v)
...
AttributeError: 'MutableDenseMatrix' object has no attribute 'is_commutative'
>>> outer(A*u, v)
...
AttributeError: 'MutableDenseMatrix' object has no attribute 'is_commutative'

>>> D = Matrix([[div(u), 0, 0], [0, div(u), 0], [0, 0, div(u)]])

>>> inner(D,D)
...
AttributeError: 'MutableDenseMatrix' object has no attribute 'is_commutative'

The problem comes from the same piece of code when treating the Mul objects;

        if isinstance(arg1, Mul):
            a = arg1.args
        else:
            a = [arg1]

        if isinstance(arg2, Mul):
            b = arg2.args
        else:
            b = [arg2]

        args_1 = [i for i in a if not i.is_commutative]
        c1     = [i for i in a if not i in args_1]
        args_2 = [i for i in b if not i.is_commutative]
        c2     = [i for i in b if not i in args_2]

        a = reduce(mul, args_1)
        b = reduce(mul, args_2)
        c = Mul(*c1)*Mul(*c2)

[e-moral-sanchez]

We want to be able to assemble the bilinear form $\int_{\Omega} (\boldsymbol{b} \times \boldsymbol{u}) \cdot \boldsymbol{v} d \boldsymbol{x}$ using the cross product matrix.
Fixing this bug will allow us to use the matrix as intended BilinearForm((u, v), integral(domain, dot(S*u, v))), instead of having to type the expression component by component.