Devel matrices

sympde/calculus/__init__.py

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

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,13 @@
 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
+
+
 class MatrixSymbolicExpr(Expr):
     is_commutative = False
     _op_priority   = 20.0
@@ -289,6 +293,124 @@ def _sympystr(self, printer):
         sstr = printer.doprint
         return '{}[{}]'.format(sstr(self.args[0]),sstr(self.args[1]))
 
+class Matrix(MatrixSymbolicExpr):
+
+    def __new__(cls, *args, **options):
+        name = options.pop('name')
+        if name is None:
+            raise ValueError('Expecting a name keyword.')
+
+        if not( len(args) == 1 ):
+            raise ValueError('Expecting one argument.')
+
+        _args = args[0]
+        if not isinstance(_args, list):
+            raise TypeError('Expecting a list.')
+
+        for _ in _args:
+            if not isinstance(_, list):
+                raise TypeError('args components must be a list.')
+
+        sizes = [len(_) for _ in _args]
+        if not( len(set(sizes)) == 1):
+            raise ValueError('Wrong matrix size')
+
+        # make _args a Tuple of Tuples
+        newargs = []
+        for _ in _args:
+            a = Tuple(*_)
+            newargs.append(a)
+        _args = Tuple(*newargs)
+
+        args = (_args, args[1:])
+
+        obj = Expr.__new__(cls, *args)
+        obj._name = name
+
+        return obj
+
+    @property
+    def name(self):
+        return self._name
+
+    def __getitem__(self, key):
+        try:
+            i,j = key
+            return self.args[0][i][j]
+        except:
+            return MatrixElement(self, key)
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return '{}'.format(sstr(self.name))
+
+    def __hash__(self):
+        return hash((self.name, self.args))
+
+    def __eq__(self, a):
+        if isinstance(a, Matrix):
+            eq = self.name == a.name
+            return eq
+        return False
+
+    def to_sympy(self):
+        return ImmutableDenseMatrix(self.args[0])
+
+
+class Vector(MatrixSymbolicExpr):
+
+    def __new__(cls, *args, **options):
+        name = options.pop('name')
+        if name is None:
+            raise ValueError('Expecting a name keyword.')
+
+        if not( len(args) == 1 ):
+            raise ValueError('Expecting one argument.')
+
+        _args = args[0]
+        if not isinstance(_args, list):
+            raise TypeError('Expecting a list.')
+
+        # make _args a Tuple
+        _args = Tuple(*_args)
+
+        args = (_args, args[1:])
+
+        obj = Expr.__new__(cls, *args)
+        obj._name = name
+
+        return obj
+
+    @property
+    def name(self):
+        return self._name
+
+    def __getitem__(self, key):
+        try:
+            i = key
+            return self.args[0][i]
+        except:
+            # TODO ARA do we keep returning a MatrixElement?
+            return MatrixElement(self, key)
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return '{}'.format(sstr(self.name))
+
+    def __hash__(self):
+        return hash((self.name, self.args))
+
+    def __eq__(self, a):
+        if isinstance(a, Matrix):
+            eq = self.name == a.name
+            return eq
+        return False
+
+    def to_sympy(self):
+        _args = [[_] for _ 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,13 @@ 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 ImmutableDenseMatrix([expr.args[0]])
+            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))