Add InterfaceMapping (#69)

InterfaceMapping is used to represent a mapping in the interface.
It's needed in LogicalExpr which takes the mapping as argument in order
to translate an expression from the physical domain to the logical domain.
So it needs to use the right mapping for a given atom in an interface expression.

The changes are :

* Add InterfaceMapping class
* Use OrderedDict in TerminalExpr instead of Dict
* Add AffineMapping class
* Add the __add__ magic method in Integral classes

sympde/expr/evaluation.py

@@ -74,7 +74,7 @@ def zero_matrix(n_rows, n_cols):
 #==============================================================================
 def _get_size_and_starts(ls):
     n = 0
-    d_indices = {}
+    d_indices = OrderedDict()
     for x in ls:
         d_indices[x] = n
         if isinstance(x, ScalarTestFunction):
@@ -260,7 +260,7 @@ def _split_expr_over_interface(expr, interface, tests=None, trials=None):
     # ...
 
     int_expressions = []
-    bnd_expressions = {}
+    bnd_expressions = OrderedDict()
 
     # ...
     # we replace all jumps
@@ -272,7 +272,7 @@ def _split_expr_over_interface(expr, interface, tests=None, trials=None):
     # ...
 
     # ...
-    d_trials = {}
+    d_trials = OrderedDict()
     for u in trials:
         u_minus = minus(u)
         u_plus  = plus(u)
@@ -281,7 +281,7 @@ def _split_expr_over_interface(expr, interface, tests=None, trials=None):
 #        # TODO add sub for avg
 #        expr = expr.subs({jump(u): u_minus - u_plus})
 
-    d_tests  = {}
+    d_tests  = OrderedDict()
     for v in tests:
         v_minus = minus(v)
         v_plus  = plus(v)
@@ -496,14 +496,14 @@ def eval(cls, *_args, **kwargs):
             dim = expr.ldim
             domain = expr.domain
             if isinstance(domain, Union):
-                domain = list(domain._args)
+                domain = domain.as_tuple()
 
             elif not is_sequence(domain):
-                domain = [domain]
+                domain = (domain,)
             # ...
 
             # ...
-            d_expr = {}
+            d_expr = OrderedDict()
             for d in domain:
                 d_expr[d] = S.Zero
             # ...
@@ -528,7 +528,6 @@ def eval(cls, *_args, **kwargs):
                     for d in domain:
                         d_expr[d] += newexpr
                     # ...
-
             else:
                 newexpr = cls.eval(expr.expr, dim=dim)
                 newexpr = expand(newexpr)
@@ -554,7 +553,7 @@ def eval(cls, *_args, **kwargs):
             # ...
 
             # ...
-            d_new = {}
+            d_new = OrderedDict()
             for domain, newexpr in d_expr.items():
 
                 if newexpr != 0:
@@ -571,7 +570,7 @@ def eval(cls, *_args, **kwargs):
 
             # ...
             ls = []
-            d_all = {}
+            d_all = OrderedDict()
 
             # ... treating interfaces
             keys = [k for k in d_new.keys() if isinstance(k, Interface)]
@@ -609,7 +608,7 @@ def eval(cls, *_args, **kwargs):
             for domain in keys:
 
                 newexpr = d_new[domain]
-                d       = {interior : newexpr for interior in domain.as_tuple()}
+                d       = OrderedDict((interior, newexpr) for interior in domain.as_tuple())
                             # ...
                 for k, v in d.items():
                     if k in d_all.keys():
@@ -618,7 +617,7 @@ def eval(cls, *_args, **kwargs):
                     else:
                         d_all[k] = v
 
-            d = {}
+            d = OrderedDict()
 
             for k, v in d_new.items():
                 if not isinstance( k, (Interface, Union)):
@@ -690,7 +689,6 @@ def eval(cls, *_args, **kwargs):
                     for i in range(0, dim):
                         e += M[i] * n[i]
                     return e
-
             else:
                 raise ValueError('> Only traces of order 0 and 1 are available')
 
@@ -955,7 +953,7 @@ def eval(cls, *_args, **kwargs):
             raise ValueError('Expecting one argument')
 
         expr = _args[0]
-        d_atoms = kwargs.pop('d_atoms', {})
+        d_atoms = kwargs.pop('d_atoms', OrderedDict())
         mapping = kwargs.pop('mapping', None)
 
         if isinstance(expr, Add):
@@ -1032,7 +1030,7 @@ def eval(cls, *_args, **kwargs):
             variables  = list(terminal_expr.atoms(ScalarTestFunction))
             variables += list(terminal_expr.atoms(IndexedTestTrial))
 
-            d_atoms = {}
+            d_atoms = OrderedDict()
             for a in variables:
                 new = _split_test_function(a)
                 d_atoms[a] = new[a]

sympde/expr/expr.py

@@ -54,16 +54,17 @@ def _get_domain(expr):
     if isinstance(expr, (DomainIntegral, BoundaryIntegral, InterfaceIntegral)):
         return expr.domain
 
-    elif isinstance(expr, (Add,Mul)):
-        domains = set()
+    elif isinstance(expr, (Add, Mul)):
+        domains = []
         for a in expr.args:
             a = _get_domain(a)
             if isinstance(a, Union):
-                domains = domains.union(a.args)
+                domains.extend(list(a.args))
             elif isinstance(a, BasicDomain):
-                domains = domains.union([a])
+                domains.append(a)
         if len(domains) == 1:
-            return tuple(domains)[0]
+            return domains[0]
+
         return Union(*domains)
 
 #==============================================================================
@@ -104,6 +105,8 @@ class Integral(CalculusFunction):
 
     def __new__(cls, expr, domain, **options):
         # (Try to) sympify args first
+        if domain is None:
+            return sy_Zero()
 
         assert isinstance(domain, BasicDomain)
         return Integral.eval(expr, domain)
@@ -117,7 +120,7 @@ def eval(expr, domain):
             raise TypeError('only Expr are accepted')
 
         if isinstance(expr, sy_Zero):
-            return sy_Zero
+            return sy_Zero()
 
         if isinstance(expr, Add):
             args = [Integral.eval(a, domain) for a in expr.args]
@@ -152,6 +155,18 @@ def __mul__(self, o):
     def __rmul__(self, o):
         return DomainIntegral(self.expr*o, self.domain)
 
+    def __add__(self, o):
+        if o == 0:
+            return self
+        else:
+            return Add(self, o)
+
+    def __radd__(self, o):
+        if o == 0:
+            return self
+        else:
+            return Add(self, o)
+
     def __eq__(self, a):
         if isinstance(a, DomainIntegral):
             eq = self.domain == a.domain
@@ -216,6 +231,18 @@ def __mul__(self, o):
     def __rmul__(self, o):
         return BoundaryIntegral(self.expr*o, self.domain)
 
+    def __add__(self, o):
+        if o == 0:
+            return self
+        else:
+            return Add(self, o)
+
+    def __radd__(self, o):
+        if o == 0:
+            return self
+        else:
+            return Add(self, o)
+
     def __eq__(self, a):
         if isinstance(a, BoundaryIntegral):
             eq = self.domain == a.domain
@@ -244,6 +271,18 @@ def __mul__(self, o):
     def __rmul__(self, o):
         return InterfaceIntegral(self.expr*o, self.domain)
 
+    def __add__(self, o):
+        if o == 0:
+            return self
+        else:
+            return Add(self, o)
+
+    def __radd__(self, o):
+        if o == 0:
+            return self
+        else:
+            return Add(self, o)
+
     def __eq__(self, a):
         if isinstance(a, InterfaceIntegral):
             eq = self.domain == a.domain
@@ -309,7 +348,7 @@ def __new__(cls, arguments, expr):
 
         # Trivial case: null expression
         if expr == 0:
-            return sy_Zero
+            return sy_Zero()
 
         # Check that integral expression is given
         integral_types = DomainIntegral, BoundaryIntegral, InterfaceIntegral
@@ -401,7 +440,7 @@ def __new__(cls, arguments, expr):
 
         # Trivial case: null expression
         if expr == 0:
-            return sy_Zero
+            return sy_Zero()
 
         # Check that integral expression is given
         integral_types = DomainIntegral, BoundaryIntegral, InterfaceIntegral
@@ -410,6 +449,7 @@ def __new__(cls, arguments, expr):
 
         # Expand integral expression and sanitize arguments
         # TODO: why do we 'sanitize' here?
+
         expr = expand(expr)
         args = _sanitize_arguments(arguments, is_bilinear=True)
 

sympde/topology/basic.py

@@ -113,8 +113,8 @@ def __new__(cls, *args):
         for union in unions:
             args += list(union.as_tuple())
 
-        # Sort domains by name
-        args = sorted(args, key=lambda x: x.name)
+        # remove duplicates and Sort domains by name
+        args = sorted(set(args), key=lambda x: x.name)
 
         # a. If the required Union contains no domains, return None;
         # b. If it contains a single domain, return the domain itself;

sympde/topology/domain.py

@@ -57,10 +57,14 @@ def __new__(cls, name, interiors=None, boundaries=None, dim=None,
                 interiors = [interiors]
 
             else:
-                unions    = [i for i in interiors if isinstance(i, Union)]
-                interiors = [i for i in interiors if not isinstance(i, Union)]
-                for union in unions:
-                    interiors += list(union.as_tuple())
+                new_interiors = []
+                for i in interiors:
+                    if isinstance(i , Union):
+                        new_interiors += list(i.as_tuple())
+                    else:
+                        new_interiors.append(i)
+
+                interiors = new_interiors
 
                 if not all([isinstance(i, InteriorDomain) for i in interiors]):
                     raise TypeError('> all interiors must be of type InteriorDomain')
@@ -345,15 +349,16 @@ def from_file( cls, filename ):
                               boundaries=boundary,
                               connectivity=connectivity)
 
-    def join(self, other, name, bnd_minus, bnd_plus):
+    def join(self, other, name, bnd_minus=None, bnd_plus=None):
         # ... interiors
         interiors = [self.interior, other.interior]
         # ...
 
         # ... connectivity
         connectivity = Connectivity()
         # TODO be careful with '|' in psydac
-        connectivity['{l}|{r}'.format(l=bnd_minus.domain.name, r=bnd_plus.domain.name)] = (bnd_minus, bnd_plus)
+        if bnd_minus and bnd_plus:
+            connectivity['{l}|{r}'.format(l=bnd_minus.domain.name, r=bnd_plus.domain.name)] = (bnd_minus, bnd_plus)
         for k,v in self.connectivity.items():
             connectivity[k] = v