Devel matrices

sympde/calculus/__init__.py

@@ -1,3 +1,2 @@
 from .core     import *
 from .errors   import *
-from .matrices import *

sympde/calculus/core.py

@@ -428,6 +428,14 @@ def __new__(cls, arg1, arg2, **options):
             a,b = b,a
             c   = -c
 
+        # after simplifications, the expression can be reduced to 0 in some cases
+        # here's an example;
+        # if b is a Vector and u and v are scalar functions,
+        # we should have the following
+        # cross(b*u,b*v) = u*v*cross(b,b) = 0
+        if a == b:
+            return S.Zero
+
         obj = Basic.__new__(cls, a, b)
 
         return c*obj

sympde/core/__init__.py

@@ -1,3 +1,4 @@
 from .algebra      import *
 from .basic        import *
 from .utils        import *
+from .matrices     import *

sympde/calculus/matrices.py → sympde/core/matrices.py

@@ -1,9 +1,15 @@
 from sympy                 import Expr, S
 from sympy                 import Add, Mul, Pow
+from sympy                 import Tuple
 from sympy                 import sympify
 from sympy.core.decorators import call_highest_priority
 from sympde.core.basic     import _coeffs_registery, Basic
 
+from sympy import ImmutableDenseMatrix
+
+import functools
+
+
 class MatrixSymbolicExpr(Expr):
     is_commutative = False
     _op_priority   = 20.0
@@ -289,6 +295,106 @@ def _sympystr(self, printer):
         sstr = printer.doprint
         return '{}[{}]'.format(sstr(self.args[0]),sstr(self.args[1]))
 
+class Matrix(MatrixSymbolicExpr):
+    def __new__(cls, mat, *, name):
+        if not isinstance(mat, list):
+            raise TypeError('Positional argument `mat` should be a list of lists.')
+
+        for row in mat:
+            if not isinstance(row, list):
+                raise TypeError('Each row of `mat` should be a list.')
+
+        row_sizes = {len(row) for row in mat}
+        if len(row_sizes) != 1:
+            raise ValueError('Each row of `mat` should have the same length.')
+
+        if not isinstance(name, str):
+            raise TypeError('Keyword-only argument `name` must be a string.')
+
+        # Call constructor of base class Expr with immutable arguments
+        new_mat = Tuple(*[Tuple(*row) for row in mat])
+        obj = Expr.__new__(cls, new_mat, name)
+        return obj
+
+    @property
+    def name(self):
+        return self.args[1]
+
+    def __getitem__(self, key):
+        i,j = key
+        return self.args[0][i][j]
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return '{}'.format(sstr(self.name))
+
+    def __hash__(self):
+        return hash((self.name, self.args))
+
+    def __eq__(self, other):
+        if not isinstance(other, Matrix):
+            return False
+
+        if not self.name == other.name:
+            return False
+
+        result = all(x == y for x, y in zip(self.args[0], other.args[0]))
+        return result
+
+    def to_sympy(self):
+        return ImmutableDenseMatrix(self.args[0])
+
+
+class Vector(MatrixSymbolicExpr):
+
+    def __new__(cls, vec, *, name):
+        if not isinstance(vec, list):
+            raise TypeError('Positional argument `vec` should be a list.')
+
+        if not isinstance(name, str):
+            raise TypeError('Keyword-only argument `name` must be a string.')
+
+        # Call constructor of base class Expr with immutable arguments
+        new_vec = Tuple(*vec)
+        obj = Expr.__new__(cls, new_vec, name)
+        return obj
+
+    @property
+    def name(self):
+        return self.args[1]
+
+    def __getitem__(self, key):
+        i = key
+        return self.args[0][i]
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return '{}'.format(sstr(self.name))
+
+    def __hash__(self):
+        return hash((self.name, self.args))
+
+    def __eq__(self, other):
+        # TODO BUG
+        # We should check that other is a Vector
+        # right now, there is a problem when using Cross
+        # see the linearity test of
+        # f = lambda u,v: cross(b*u, b*v)
+        # where b is a Vector
+        if not isinstance(other, Vector):
+            return False
+
+        if not self.name == other.name:
+            return False
+
+        result = self.args[0] == other.args[0]
+        return result
+
+    def to_sympy(self):
+        args = [[a] for a in self.args[0]]
+        return ImmutableDenseMatrix(args)
+
+
 Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = {
     "Mul": [lambda x: MatSymbolicMul(*x.args)],
     "Add": [lambda x: MatSymbolicAdd(*x.args)]

sympde/expr/evaluation.py

@@ -20,12 +20,14 @@
                                  Dot_2d, Inner_2d, Cross_2d,
                                  Dot_3d, Inner_3d, Cross_3d)
 from sympde.core.utils import random_string
+from sympde.core.matrices import SymbolicDeterminant, Inverse, Transpose
+from sympde.core.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
+from sympde.core.matrices import Matrix as SympdeMatrix
+from sympde.core.matrices import Vector as SympdeVector
 
 from sympde.calculus import jump, avg, minus, plus
 from sympde.calculus import Jump, is_zero
 from sympde.calculus.core import _generic_ops, _diff_ops
-from sympde.calculus.matrices import SymbolicDeterminant, Inverse, Transpose
-from sympde.calculus.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
 
 from sympde.topology.basic   import BasicDomain, Union, Interval
 from sympde.topology.basic   import Boundary, Interface
@@ -614,6 +616,12 @@ def eval(cls, expr, domain):
         elif isinstance(expr, Inverse):
             return cls.eval(expr.arg, domain=domain).inv()
 
+        elif isinstance(expr, SympdeMatrix):
+            return expr.to_sympy()
+
+        elif isinstance(expr, SympdeVector):
+            return expr.to_sympy()
+
         elif isinstance(expr, ScalarFunction):
             return expr
 

sympde/expr/expr.py

@@ -22,6 +22,8 @@
 from sympde.topology.space import ScalarFunction
 from sympde.topology.space import VectorFunction
 from sympde.topology.space import Trace, trace_0, trace_1
+from sympde.core.matrices import Matrix as SympdeMatrix
+from sympde.core.matrices import Vector
 
 from .errors import UnconsistentLinearExpressionError
 from .basic  import BasicForm
@@ -538,9 +540,12 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
         # ...
 
         # ...
-        is_vector = isinstance(expr, (Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
+        is_vector = isinstance(expr, (Vector, SympdeMatrix, Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
         if is_vector:
-            expr = ImmutableDenseMatrix(expr)
+            if isinstance(expr, (Vector, SympdeMatrix)):
+                expr = expr.to_sympy()
+            else:
+                expr = ImmutableDenseMatrix(expr)
         # ...
 
         # ...
@@ -610,9 +615,12 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
         # ...
 
         # ...
-        is_vector = isinstance(expr, (Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
+        is_vector = isinstance(expr, (Vector, SympdeMatrix, Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
         if is_vector:
-            expr = ImmutableDenseMatrix(expr)
+            if isinstance(expr, (Vector, SympdeMatrix)):
+                expr = expr.to_sympy()
+            else:
+                expr = ImmutableDenseMatrix(expr)
         # ...
 
         # ...

sympde/expr/tests/test_equation_2d.py

@@ -2,7 +2,6 @@
 
 #import pytest
 
-from sympy.core.containers import Tuple
 from sympy import pi, cos, sin
 from sympy import ImmutableDenseMatrix as Matrix
 
@@ -48,7 +47,7 @@ def test_equation_2d_1():
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
     int_2 = lambda expr: integral(B2, expr)
-    
+
     # ... bilinear/linear forms
     expr = dot(grad(v), grad(u))
     a1 = BilinearForm((v,u), int_0(expr))
@@ -169,7 +168,7 @@ def test_equation_2d_2():
     alpha = Constant('alpha', real=True)
 
     int_0 = lambda expr: integral(domain , expr)
-    
+
     s  = BilinearForm((tau,sigma), int_0(dot(grad(tau), grad(sigma))))
     m  = BilinearForm((tau,sigma), int_0(tau*sigma))
     b1 = BilinearForm((tau,dw), int_0(bracket(pn, dw) * tau))
@@ -202,7 +201,7 @@ def test_equation_2d_3():
     x,y = domain.coordinates
 
     B1 = Boundary(r'\Gamma_1', domain)
-    
+
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
@@ -235,7 +234,7 @@ def test_equation_2d_4():
     x,y = domain.coordinates
 
     B1 = Boundary(r'\Gamma_1', domain)
-    
+
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
@@ -283,7 +282,7 @@ def test_equation_2d_5():
     p,q = [element_of(V, name=i) for i in ['p', 'q']]
 
     int_0 = lambda expr: integral(domain , expr)
-    
+
     a0 = BilinearForm((v,u), int_0(inner(grad(v), grad(u))))
     print('     a0 done.')
     a1 = BilinearForm((q,p), int_0(p*q))
@@ -335,7 +334,7 @@ def test_equation_2d_6():
     u,v = [element_of(V, name=i) for i in ['u', 'v']]
 
     int_0 = lambda expr: integral(domain , expr)
-    
+
     # ...
     expr = kappa * dot(grad(u), grad(v)) + dot(b, grad(u)) * v
     a = BilinearForm((v,u), int_0(expr))