Some bug-fixes for Psydac (#91)

In order to avoid problems when assembling vector expressions over the boundary, we do not create `Trace` objects from the expressions. Such a kind of reverse engineering is complex and error prone, and not really needed at the moment.

* Fix bug in Integral.subs_boundary_expr
* Set Trace commutativity from its expression
* Let IndexedTestTrial use indices of type sympy.Integer
* Don't substitute u with Trace(u) in boundary integral
* Deactivate unit tests for Trace
* Update version.py

sympde/expr/expr.py

@@ -154,12 +154,19 @@ def __new__(cls, expr, domain, **options):
         if isinstance(domain, Union):
             exprs = [cls(expr, d) for d in domain]
             return IntAdd(*exprs)
+
         elif isinstance(domain, Domain):
             interiors = domain.interior if isinstance(domain.interior, Union) else [domain.interior]
             exprs = [cls(expr, d) for d in interiors]
             return IntAdd(*exprs)
+
         elif isinstance(domain, Boundary):
-            expr = cls.subs_boundary_expr(expr, domain)
+            #------------------------------------------------------------
+            # NOTE [YG, 27.01.2021]:
+            #   we stop using Traces in the boundary integrals because of
+            #   an error that arises when using cross(v, nn).
+            #------------------------------------------------------------
+#            expr = cls.subs_boundary_expr(expr, domain)
             obj = CalculusFunction.__new__(cls, expr, domain)
             obj.is_boundary_integral = True
 
@@ -170,6 +177,7 @@ def __new__(cls, expr, domain, **options):
         elif isinstance(domain, InteriorDomain):
             obj = CalculusFunction.__new__(cls, expr, domain)
             obj.is_domain_integral = True
+
         else:
             raise TypeError(domain)
 
@@ -228,14 +236,10 @@ def __hash__(self):
 
     @classmethod
     def subs_boundary_expr(cls, expr, domain):
-        atoms_1 = list(expr.atoms(Dot,Trace))
-
-        for i in range(len(atoms_1)):
-            a = atoms_1[i]
-            if isinstance(a, Dot):
-                if not isinstance(a.args[0], NormalVector):
-                    if not isinstance(a.args[1], NormalVector):
-                        atoms_1.remove(a)
+
+        atoms_1 = [*expr.atoms(Trace)] + [a for a in expr.atoms(Dot) \
+                                       if isinstance(a.args[0], NormalVector) \
+                                       or isinstance(a.args[1], NormalVector)]
 
         subs_1  = {a:Dummy() for a in atoms_1}
         expr, _ = expr._xreplace(subs_1)

sympde/expr/tests/test_expr.py

@@ -22,7 +22,7 @@
 from sympde.topology import InteriorDomain, Union
 from sympde.topology import Boundary, NormalVector
 from sympde.topology import Domain
-from sympde.topology import trace_1
+#from sympde.topology import trace_1  # TODO [YG, 27.01.2021]: fix trace
 from sympde.topology import Square
 from sympde.topology import ElementDomain
 from sympde.topology import Area
@@ -137,7 +137,8 @@ def test_linear_form_2d_1():
     print(l)
 
     assert(l.domain == B1)
-    assert(l(v1) == int_1(v1*trace_1(g, B1)))
+#    assert(l(v1) == int_1(v1*trace_1(g, B1))) # TODO [YG, 27.01.2021]: fix trace
+    assert(l(v1) == int_1(v1*dot(g, nn)))
     # ...
 
     # ...
@@ -148,7 +149,8 @@ def test_linear_form_2d_1():
     for i in l.domain.args:
         assert(i in [domain.interior, B1])
 
-    assert(l(v1) == int_1(v1*trace_1(g, B1)) + int_0(x*y*v1))
+#    assert(l(v1) == int_1(v1*trace_1(g, B1)) + int_0(x*y*v1)) # TODO [YG, 27.01.2021]: fix trace
+    assert(l(v1) == int_1(v1*dot(g, nn)) + int_0(x*y*v1))
     # ...
 
     # ...
@@ -277,7 +279,8 @@ def test_bilinear_form_2d_1():
     a = BilinearForm((u,v), int_1(v*dot(grad(u), nn)))
 
     assert(a.domain == B1)
-    assert(a(u1,v1) == int_1(v1*trace_1(grad(u1), B1)))
+#    assert(a(u1,v1) == int_1(v1*trace_1(grad(u1), B1))) # TODO [YG, 27.01.2021]: fix trace
+    assert a(u1,v1) == int_1(v1*dot(grad(u1), nn))
     # ...
 
     # ...
@@ -288,7 +291,8 @@ def test_bilinear_form_2d_1():
     for i in a.domain.args:
         assert(i in [domain.interior, B1])
 
-    assert(a(u1,v1) == int_0(u1*v1) + int_1(v1*trace_1(grad(u1), B1)))
+#    assert(a(u1,v1) == int_0(u1*v1) + int_1(v1*trace_1(grad(u1), B1))) # TODO [YG, 27.01.2021]: fix trace
+    assert a(u1,v1) == int_0(u1*v1) + int_1(v1*dot(grad(u1), nn))
     # ...
 
     # ...

sympde/topology/space.py

@@ -516,10 +516,9 @@ def __new__(cls, base, *args, **kw_args):
     @property
     def free_symbols(self):
         base_free_symbols = self.base.free_symbols
-        symbolic_indices = [i for i in self.indices if isinstance(i, Basic)]
-        if len(symbolic_indices) > 0:
-            raise ValueError('symbolic indices not yet available')
-
+        for i in self.indices:
+            if not isinstance(i, (int, Integer)):
+                raise ValueError('Symbolic index {} of type {} cannot be used'.format(i, type(i)))
         return base_free_symbols
 
     @property
@@ -778,18 +777,22 @@ class Trace(AtomicExpr):
     Represents the trace over a boundary and a space function
 
     """
-    is_commutative = True
+    is_commutative = None
+
     def __new__(cls, expr, boundary, order=0, **options):
 #        # TODO these tests are not working for the moment for Grad(u)
 #        if not expr.atoms((ScalarTestFunction, VectorTestFunction, ScalarField)):
 #            raise TypeError('> Wrong type for expr')
 #
 #        if not(expr.space.domain is boundary.domain):
 #            raise ValueError('> Space and boundary domains must be the same')
-        evaluate = options.pop('evaluate',True)
-        if evaluate:
+
+        if options.pop('evaluate',True):
             return cls.eval(expr, boundary, order)
-        return Basic.__new__(cls, expr, boundary, order)
+
+        obj = Basic.__new__(cls, expr, boundary, order)
+        obj.is_commutative = expr.is_commutative
+        return obj
 
     @property
     def expr(self):

sympde/version.py

@@ -1 +1 @@
-__version__ = "0.10.7"
+__version__ = "0.10.8"