Devel mapping (#82)

The following changes have been made:
- Add computation of Jacobian matrix and its inverse in creation of mapping;
- Jacobian matrix and its inverse may be explicitly provided by user when creating symbolic mapping;
- Add CallableMapping where we make symbolic mapping callable;
- Allow for pdim (number of physical dimensions) greater than ldim (number of logical dimensions).

sympde/calculus/core.py

@@ -125,6 +125,7 @@ def has(obj, types):
 
 @cacheit
 def is_zero(x):
+
     if isinstance(x, (Matrix, ImmutableDenseMatrix)):
         return all( i==0 for i in x[:])
     else:

sympde/calculus/matrices.py

@@ -1,14 +1,16 @@
 from sympy                 import Expr, S
 from sympy                 import Add, Mul, Pow
-from sympy.core.decorators import call_highest_priority
-from sympy                 import Basic
 from sympy                 import sympify
-from sympde.core.basic     import _coeffs_registery
+from sympy.core.decorators import call_highest_priority
+from sympde.core.basic     import _coeffs_registery, Basic
 
 class MatrixSymbolicExpr(Expr):
     is_commutative = False
     _op_priority   = 20.0
     is_MatrixSymbolicExpr = True
+    is_Matrix     = True
+    is_ZeroMatrix = False
+    is_Identity   = False
 
     def __new__(cls, *args, **options):
         return Expr.__new__(cls, *args)
@@ -34,7 +36,7 @@ def trace(self):
 
     # The following is adapted from the core Expr object
     def __neg__(self):
-        return MatMul(S.NegativeOne, self)
+        return MatSymbolicMul(S.NegativeOne, self)
 
     def __abs__(self):
         return MatAbs(self)
@@ -44,39 +46,39 @@ def __getitem__(self, key):
 
     @call_highest_priority('__radd__')
     def __add__(self, other):
-        return MatAdd(self, other)
+        return MatSymbolicAdd(self, other)
 
     @call_highest_priority('__add__')
     def __radd__(self, other):
-        return MatAdd(other, self)
+        return MatSymbolicAdd(other, self)
 
     @call_highest_priority('__rsub__')
     def __sub__(self, other):
-        return MatAdd(self, -other)
+        return MatSymbolicAdd(self, -other)
 
     @call_highest_priority('__sub__')
     def __rsub__(self, other):
-        return MatAdd(other, -self)
+        return MatSymbolicAdd(other, -self)
 
     @call_highest_priority('__rmul__')
     def __mul__(self, other):
-        return MatMul(self, other)
+        return MatSymbolicMul(self, other)
 
     @call_highest_priority('__mul__')
     def __rmul__(self, other):
-        return MatMul(other, self)
+        return MatSymbolicMul(other, self)
 
     @call_highest_priority('__pow__')
     def __pow__(self, other):
-        return MatPow(self, other)
+        return MatSymbolicPow(self, other)
 
     @call_highest_priority('__div__')
     def __div__(self, other):
-        return MatMul(self , other**S.NegativeOne)
+        return MatSymbolicMul(self , other**S.NegativeOne)
 
     @call_highest_priority('__rdiv__')
     def __rdiv__(self, other):
-        return MatMul(other, Pow(self, S.NegativeOne))
+        return MatSymbolicMul(other, Pow(self, S.NegativeOne))
 
     __truediv__ = __div__
     __rtruediv__ = __rdiv__
@@ -113,8 +115,8 @@ def _sympystr(self, printer):
         sstr = printer.doprint
         return '{}.T'.format(sstr(self.arg))
 
-class MatMul(MatrixSymbolicExpr, Mul):
-    is_MatMul = True
+class MatSymbolicMul(MatrixSymbolicExpr, Mul):
+
     def __new__(cls, *args, **options):
         args = [sympify(a) for a in args if a != 1]
 
@@ -123,7 +125,7 @@ def __new__(cls, *args, **options):
                 newargs = []
                 for e in a.args:
                     t = args[:i] + [e] + args[i+1:]
-                    newargs.append(MatMul(*t))
+                    newargs.append(MatSymbolicMul(*t))
                 return type(a)(*newargs)
 
         newargs = []
@@ -174,19 +176,19 @@ def _expandsums(sums):
         if L == 1:
             return sums[0].args
         terms = []
-        left = MatMul._expandsums(sums[:L//2])
-        right = MatMul._expandsums(sums[L//2:])
-        terms = [MatMul(a, b) for a in left for b in right]
-        added = MatAdd(*terms)
+        left = MatSymbolicMul._expandsums(sums[:L//2])
+        right = MatSymbolicMul._expandsums(sums[L//2:])
+        terms = [MatSymbolicMul(a, b) for a in left for b in right]
+        added = MatSymbolicAdd(*terms)
         return Add.make_args(added) # it may have collapsed down to one term
 
-class MatAdd(MatrixSymbolicExpr, Add):
-    is_MatAdd = True
+class MatSymbolicAdd(MatrixSymbolicExpr, Add):
+
     def __new__(cls, *args, **options):
         args = [sympify(a) for a in args if a != 0]
         newargs = []
         for i in args:
-            if isinstance(i, MatAdd):
+            if isinstance(i, MatSymbolicAdd):
                 newargs += list(i.args)
             else:
                 newargs.append(i)
@@ -203,15 +205,15 @@ def __new__(cls, *args, **options):
     def _sympystr(self, printer):
         sstr = printer.doprint
         #return ' + '.join((sstr(i) for i in self.args))
-        return 'MatAdd({})'.format(','.join((sstr(i) for i in self.args)))
+        return 'MatSymbolicAdd({})'.format(','.join((sstr(i) for i in self.args)))
 
-class MatPow(MatrixSymbolicExpr, Pow):
+class MatSymbolicPow(MatrixSymbolicExpr, Pow):
     def _sympystr(self, printer):
         sstr = printer.doprint
         #return '**'.join((sstr(i) for i in self.args))
-        return 'MatPow({})'.format(','.join((sstr(i) for i in self.args)))
+        return 'MatSymbolicPow({})'.format(','.join((sstr(i) for i in self.args)))
 
-class MatAbs(MatrixSymbolicExpr):
+class MatSymbolicAbs(MatrixSymbolicExpr):
     def _sympystr(self, printer):
         sstr = printer.doprint
         return 'Abs({})'.format(sstr(self.args[0]))
@@ -272,3 +274,13 @@ def _sympystr(self, printer):
         sstr = printer.doprint
         return '{}[{}]'.format(sstr(self.args[0]),sstr(self.args[1]))
 
+def f1(*x):
+    return MatSymbolicMul(*x)
+def f2(*x):
+    return MatSymbolicAdd(*x)
+
+Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = {
+    "Mul": [f1],
+    "Add": [f2]
+}
+

sympde/expr/evaluation.py

@@ -23,9 +23,9 @@
 from sympde.calculus.core import _generic_ops, _diff_ops
 
 from sympde.calculus.matrices import SymbolicDeterminant, Inverse, Transpose
-from sympde.calculus.matrices import MatPow, MatrixElement, SymbolicTrace
+from sympde.calculus.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
 
-from sympde.topology.mapping import Jacobian, JacobianSymbol, InterfaceMapping
+from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, MultiPatchMapping, JacobianInverseSymbol
 
 from sympde.topology.basic   import BasicDomain, Union, Interval
 from sympde.topology.domain  import NormalVector, TangentVector
@@ -403,6 +403,7 @@ def __new__(cls, *args, **options):
             r = cls.eval(*args, **options)
         else:
             r = None
+
         if r is None:
             return Basic.__new__(cls, *args, **options)
         else:
@@ -432,7 +433,6 @@ def _annotate(*args):
                 new_fields = [f.space.field(f.name) for f in fields]
                 expr = expr.subs(zip(fields, new_fields))
 
-
         args[0] = expr
         return args
 
@@ -477,7 +477,15 @@ def eval(cls, *_args, **kwargs):
 
         elif isinstance(expr, JacobianSymbol):
             axis = expr.axis
-            J    = Jacobian(expr.mapping)
+            J    = expr.mapping.jacobian_expr
+            if axis is None:
+                return J
+            else:
+                return J.col_del(axis)
+
+        elif isinstance(expr, JacobianInverseSymbol):
+            axis = expr.axis
+            J    = expr.mapping.jacobian_inv_expr
             if axis is None:
                 return J
             else:
@@ -495,9 +503,12 @@ def eval(cls, *_args, **kwargs):
         elif isinstance(expr, Inverse):
             return cls.eval(expr.arg, dim=dim, logical=logical).inv()
 
-        elif isinstance(expr, (ScalarTestFunction, VectorTestFunction)):
+        elif isinstance(expr, ScalarTestFunction):
             return expr
 
+        elif isinstance(expr, VectorTestFunction):
+            return ImmutableDenseMatrix([[expr[i]] for i in range(dim)])
+
         elif isinstance(expr, PullBack):
             return cls.eval(expr.expr, dim=dim, logical=True)
 
@@ -520,7 +531,7 @@ def eval(cls, *_args, **kwargs):
             # ...
             if isinstance(expr.expr, Add):
                 for a in expr.expr.args:
-                    newexpr = cls.eval(a, dim=dim, logical=True)
+                    newexpr = cls.eval(a, dim=dim, logical=logical)
 
                     # ...
                     try:
@@ -646,7 +657,8 @@ def eval(cls, *_args, **kwargs):
         elif isinstance(expr, Integral):
             dim     = expr.domain.dim if dim is None else dim
             logical = expr.domain.mapping is None
-            return cls.eval(expr._args[0], dim=dim, logical=logical)
+            expr    = cls.eval(expr._args[0], dim=dim, logical=logical)
+            return expr
 
         elif isinstance(expr, NormalVector):
             lines = [[expr[i] for i in range(dim)]]
@@ -714,10 +726,8 @@ def eval(cls, *_args, **kwargs):
 
         elif isinstance(expr, LogicalExpr):
             M         = expr.mapping
-            dim       = expr.dim
-            expr      = cls(expr.expr, dim=dim)
-            dim       = M.rdim
-            return LogicalExpr(expr, mapping=M, dim=dim)
+            expr      = cls(expr.expr, dim=expr.dim)
+            return LogicalExpr(expr, mapping=M, dim=M.ldim)
         return expr
 
 

sympde/expr/expr.py

@@ -217,8 +217,11 @@ def __eq__(self, a):
             return eq
         return False
 
+    def _hashable_content(self):
+        return (self.expr , self.domain)
+
     def __hash__(self):
-        return hash(self.expr) + hash(self.domain)
+        return hash(self._hashable_content())
 
     @classmethod
     def subs_boundary_expr(cls, expr, domain):

sympde/expr/tests/test_expr.py

@@ -861,7 +861,7 @@ def test_terminal_expr_bilinear_2d_4():
 def test_terminal_expr_bilinear_3d_1():
 
     domain = Domain('Omega', dim=3)
-    M      = Mapping('M', 3)
+    M      = Mapping('M', dim=3)
 
     mapped_domain = M(domain)
 

sympde/expr/tests/test_tensor_2d.py

@@ -73,7 +73,7 @@ def test_tensorize_2d_1_mapping():
 
     DIM = 2
 
-    M = Mapping('Map', DIM)
+    M = Mapping('Map', dim=DIM)
 
     domain = M(Domain('Omega', dim=DIM))
     B1     = Boundary(r'\Gamma_1', domain)
@@ -102,7 +102,7 @@ def test_tensorize_2d_1_mapping():
 def test_tensorize_2d_2_mapping():
 
     DIM = 2
-    M = Mapping('M', DIM)
+    M = Mapping('M', dim=DIM)
     domain = M(Domain('Omega', dim=DIM))
 
     V = VectorFunctionSpace('V', domain)

sympde/topology/__init__.py

@@ -1,7 +1,9 @@
-from .basic        import *
-from .derivatives  import *
-from .datatype     import *
-from .domain       import *
-from .mapping      import *
-from .measure      import *
-from .space        import *
+from .basic              import *
+from .derivatives        import *
+from .datatype           import *
+from .domain             import *
+from .mapping            import *
+from .measure            import *
+from .space              import *
+from .analytical_mapping import *
+from .callable_mapping   import *