Multi-patch mapping (#80)

- Add MultiPatchMapping class to handle mappings in the multi-patch case
- Fix bug in Grad and LogicalGrad for the VectorTestFunction case
- Support sympy-1.6
- Update version.py

setup.py

@@ -38,7 +38,7 @@
 
 # Dependencies
 install_requires = [
-    'sympy>=1.2,<1.6',
+    'sympy>=1.2',
     'h5py',
     'pytest',
     'pyyaml',

sympde/calculus/core.py

@@ -165,7 +165,7 @@ def __hash__(self):
 
 class BasicOperatorAdd(Add):
     _op_priority   = 10.005
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
 
         newargs = []
         for i in args:

sympde/calculus/matrices.py

@@ -10,7 +10,7 @@ class MatrixSymbolicExpr(Expr):
     _op_priority   = 20.0
     is_MatrixSymbolicExpr = True
 
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         return Expr.__new__(cls, *args)
 
     def inv(self):
@@ -83,7 +83,7 @@ def __rdiv__(self, other):
 
 class Inverse(MatrixSymbolicExpr):
     is_commutative = False
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         assert len(args) == 1
         if isinstance(args[0], Inverse):
             return args[0].arg
@@ -99,7 +99,7 @@ def _sympystr(self, printer):
 
 class Transpose(MatrixSymbolicExpr):
     is_commutative = False
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         assert len(args) == 1
         if isinstance(args[0], Transpose):
             return args[0].arg
@@ -115,7 +115,7 @@ def _sympystr(self, printer):
 
 class MatMul(MatrixSymbolicExpr, Mul):
     is_MatMul = True
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         args = [sympify(a) for a in args if a != 1]
 
         for i,a in enumerate(args):
@@ -182,7 +182,7 @@ def _expandsums(sums):
 
 class MatAdd(MatrixSymbolicExpr, Add):
     is_MatAdd = True
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         args = [sympify(a) for a in args if a != 0]
         newargs = []
         for i in args:
@@ -218,7 +218,7 @@ def _sympystr(self, printer):
 
 class SymbolicDeterminant(Expr):
     is_commutative = True
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         assert len(args) == 1
         return Basic.__new__(cls, *args)
 
@@ -233,7 +233,7 @@ def _sympystr(self, printer):
 
 class SymbolicTrace(Expr):
     is_commutative = True
-    def __new__(cls, *args):
+    def __new__(cls, *args, **options):
         assert len(args) == 1
         arg = args[0]
         if isinstance(arg, Add):
@@ -257,7 +257,7 @@ def _sympystr(self, printer):
         return 'tr({})'.format(arg)
 
 class MatrixElement(Expr):
-    def __new__(cls, base, indices):
+    def __new__(cls, base, indices, **options):
         return Expr.__new__(cls, base, indices)
 
     @property

sympde/expr/expr.py

@@ -34,8 +34,8 @@
 from functools import reduce
 
 class IntAdd(Add):
-    _op_priority          = 20
-    def __new__(cls, *args):
+    _op_priority  = 20
+    def __new__(cls, *args, **options):
         newargs = []
         for i in args:
             if isinstance(i, IntAdd):
@@ -104,7 +104,7 @@ def _get_domain(expr):
 class LinearExpr(BasicExpr):
     is_linear = True
 
-    def __new__(cls, arguments, expr):
+    def __new__(cls, arguments, expr, **options):
 
         if expr.atoms(Integral):
             raise TypeError('')
@@ -264,7 +264,7 @@ def _sympystr(self, printer):
 class Functional(BasicForm):
     is_functional = True
 
-    def __new__(cls, expr, domain, evaluate=True):
+    def __new__(cls, expr, domain, evaluate=True, **options):
 
         # compute dim from fields if available
         ls = tuple(expr.atoms(ScalarField, VectorField, ScalarTestFunction, VectorTestFunction))
@@ -308,7 +308,7 @@ def space(self):
 class LinearForm(BasicForm):
     is_linear = True
 
-    def __new__(cls, arguments, expr):
+    def __new__(cls, arguments, expr, **options):
 
         # Trivial case: null expression
         if expr == 0:
@@ -388,7 +388,7 @@ class BilinearForm(BasicForm):
     is_bilinear = True
     _is_symmetric = None
 
-    def __new__(cls, arguments, expr):
+    def __new__(cls, arguments, expr, **options):
 
         # Trivial case: null expression
         if expr == 0:
@@ -494,7 +494,7 @@ def __call__(self, trials, tests, **kwargs):
 class Norm(Functional):
     is_norm = True
 
-    def __new__(cls, expr, domain, kind='l2', evaluate=True):
+    def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
 #        # ...
 #        tests = expr.atoms((ScalarTestFunction, VectorTestFunction))
 #        if tests:

sympde/topology/derivatives.py

@@ -568,20 +568,15 @@ def eval(cls, *_args):
         if not _args:
             return
 
-        u = _args[0]
-
-        if isinstance(u, Tuple):
-            n = len(u)
-            lines = []
-            for i in range(0, n):
-                line = [dx(u)[0,i], dy(u)[0,i]]
-                lines.append(line)
-
-            v = ImmutableDenseMatrix(lines)
+        u  = _args[0]
+        du = (dx(u), dy(u))
 
+        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
+            lines = [list(d[:]) for d in du]
         else:
-            v = ImmutableDenseMatrix([[dx(u)], [dy(u)]])
+            lines = [[d] for d in du]
 
+        v = ImmutableDenseMatrix(lines)
         return v
 
 class Grad_3d(GradBasic):
@@ -593,19 +588,15 @@ def eval(cls, *_args):
         if not _args:
             return
 
-        u = _args[0]
-
-        if isinstance(u, Tuple):
-            n = len(u)
-            lines = []
-            for i in range(0, n):
-                line = [dx(u)[0,i], dy(u)[0,i], dz(u)[0,i]]
-                lines.append(line)
-
-            v = ImmutableDenseMatrix(lines)
+        u  = _args[0]
+        du = (dx(u), dy(u), dz(u))
 
+        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
+            lines = [list(d[:]) for d in du]
         else:
-            v = ImmutableDenseMatrix([[dx(u)], [dy(u)], [dz(u)]])
+            lines = [[d] for d in du]
+
+        v = ImmutableDenseMatrix(lines)
 
         return v
 
@@ -951,19 +942,15 @@ def eval(cls, *_args):
         if not _args:
             return
 
-        u = _args[0]
-
-        if isinstance(u, (Tuple, Matrix, ImmutableDenseMatrix)):
-            n = len(u)
-            lines = []
-            for i in range(0, n):
-                line = [dx1(u)[0,i], dx2(u)[0,i]]
-                lines.append(line)
-            v = ImmutableDenseMatrix(lines)
+        u  = _args[0]
+        du = (dx1(u), dx2(u))
 
+        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
+            lines = [list(d[:]) for d in du]
         else:
-            v = ImmutableDenseMatrix([[dx1(u)], [dx2(u)]])
+            lines = [[d] for d in du]
 
+        v = ImmutableDenseMatrix(lines)
         return v
 
 class LogicalGrad_3d(GradBasic):
@@ -974,20 +961,15 @@ def eval(cls, *_args):
         if not _args:
             return
 
-        u = _args[0]
-
-        if isinstance(u, (Tuple, Matrix, ImmutableDenseMatrix)):
-            n = len(u)
-            lines = []
-            for i in range(0, n):
-                line = [dx1(u)[0,i], dx2(u)[0,i], dx3(u)[0,i]]
-                lines.append(line)
-
-            v = ImmutableDenseMatrix(lines)
+        u  = _args[0]
+        du = (dx1(u), dx2(u), dx3(u))
 
+        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
+            lines = [list(d[:]) for d in du]
         else:
-            v = ImmutableDenseMatrix([[dx1(u)], [dx2(u)], [dx3(u)]])
+            lines = [[d] for d in du]
 
+        v = ImmutableDenseMatrix(lines)
         return v
 
 #==============================================================================
@@ -1162,7 +1144,7 @@ def eval(cls, *_args):
             raise NotImplementedError('TODO')
 
         return ImmutableDenseMatrix([[dx1(dx1(u)), dx1(dx2(u))],
-                       [dx1(dx2(u)), dx2(dx2(u))]])
+                                     [dx1(dx2(u)), dx2(dx2(u))]])
 
 class LogicalHessian_3d(HessianBasic):
 
@@ -1177,8 +1159,8 @@ def eval(cls, *_args):
             raise NotImplementedError('TODO')
 
         return ImmutableDenseMatrix([[dx1(dx1(u)), dx1(dx2(u)), dx1(dx3(u))],
-                       [dx1(dx2(u)), dx2(dx2(u)), dx2(dx3(u))],
-                       [dx1(dx3(u)), dx2(dx3(u)), dx3(dx3(u))]])
+                                     [dx1(dx2(u)), dx2(dx2(u)), dx2(dx3(u))],
+                                     [dx1(dx3(u)), dx2(dx3(u)), dx3(dx3(u))]])
 
 #==============================================================================
 

sympde/topology/domain.py

@@ -293,7 +293,7 @@ def from_file( cls, filename ):
         if not(ext == '.h5'):
             raise ValueError('> Only h5 files are supported')
         # ...
-        from sympde.topology.mapping import Mapping
+        from sympde.topology.mapping import Mapping, MultiPatchMapping
 
         h5  = h5py.File( filename, mode='r' )
         yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader )
@@ -363,7 +363,7 @@ def from_file( cls, filename ):
 
     def join(self, other, name, bnd_minus=None, bnd_plus=None):
 
-        from sympde.topology.mapping import InterfaceMapping
+        from sympde.topology.mapping import InterfaceMapping, MultiPatchMapping
         # ... connectivity
         connectivity = Connectivity()
         # TODO be careful with '|' in psydac
@@ -405,7 +405,7 @@ def join(self, other, name, bnd_minus=None, bnd_plus=None):
             for k,v in connectivity.items():
                 logical_connectivity[name] = v.logical_domain
 
-            mapping        = tuple(e.mapping for e in interiors)
+            mapping        = MultiPatchMapping({e.logical_domain: e.mapping for e in interiors})
             logical_domain = Domain(name,
                             interiors=logical_interiors,
                             boundaries=logical_boundaries,

sympde/topology/mapping.py

@@ -283,8 +283,9 @@ class MultiPatchMapping(Mapping):
 
     def __new__(cls, dic):
         assert isinstance( dic, dict)
-        return Basic.__new__(cls, dict)
+        return Basic.__new__(cls, dic)
 
+    @property
     def mappings(self):
         return self.args[0]
 
@@ -294,7 +295,21 @@ def is_analytical(self):
 
     @property
     def rdim(self):
-        return self.args[0].rdim
+        return list(self.mappings.values())[0].rdim
+
+    @property
+    def is_analytical(self):
+        return all(e.is_analytical for e in self.mappings.values())
+
+    def _eval_subs(self, old, new):
+        return self
+
+    def _eval_simplify(self, **kwargs):
+        return self
+
+    def __hash__(self):
+        return hash((*self.mappings.values(), *self.mappings.keys()))
+
 #==============================================================================
 class IdentityMapping(Mapping):
     """
@@ -528,7 +543,9 @@ def test(self):
 #==============================================================================
 class Jacobian(MappingApplication):
     """
-    This class calculates the Jacobian of a mapping
+    This class calculates the Jacobian of a mapping F
+    where [J_{F}]_{i,j} =  \frac{\partial F_{i}}{\partial x_{j}}
+    or simplify J_{F} =  (\nabla F)^T
 
     """
 
@@ -566,7 +583,7 @@ def eval(cls, F):
         elif rdim == 3:
             expr = LogicalGrad_3d(F)
 
-        return expr
+        return expr.T
 
 #==============================================================================
 class Covariant(MappingApplication):
@@ -847,7 +864,7 @@ def eval(cls, expr, mapping=None, dim=None, **options):
         elif isinstance(expr, laplace):
             arg = expr.args[0]
             v   = cls.eval(grad(arg), mapping=mapping, dim=dim)
-            v   = mapping.jacobian.inv()*grad(v)
+            v   = mapping.jacobian.inv().T*grad(v)
             return SymbolicTrace(v)
 
 #        elif isinstance(expr, hessian):

sympde/topology/tests/test_mapping.py

@@ -106,7 +106,8 @@ def test_mapping_3d():
     expected = Matrix(expected)
     diff     = cov-expected
     diff.simplify()
-    assert(diff.is_zero)
+
+    assert(diff.dot(diff).is_zero)
     # ...
 
     # ...

sympde/version.py

@@ -1 +1 @@
-__version__ = "0.10.0"
+__version__ = "0.10.1"