Feec dirichlet (#87)

The following changes were made:
  - Delete node annotation;
  - Compute boundary jacobian when dim = 1;
  - Update library version.

sympde/expr/basic.py

@@ -64,7 +64,6 @@ class BasicExpr(Expr):
     is_linear     = False
     is_bilinear   = False
     is_functional = False
-    is_annotated  = False
 
     @property
     def fields(self):
@@ -85,60 +84,13 @@ def constants(self):
         # no redanduncy
         return tuple(ls)
 
-    def annotate(self):
-
-        if self.is_annotated:
-            return self
-
-        if self.is_bilinear or self.is_linear:
-
-            fields          = self.fields
-            scalar_fields   = [f for f in fields if isinstance(f.space, ScalarFunctionSpace)]
-            vector_fields   = [f for f in fields if isinstance(f.space, VectorFunctionSpace)]
-
-            new_scalar_fields   = [f.space.field(f.name) for f in scalar_fields]
-            new_vector_fields   = [f.space.field(f.name) for f in vector_fields]
-
-            indexed             = self.expr.atoms(IndexedTestTrial)
-            indexed             = [f for f in indexed if f.base in vector_fields]
-            new_indexed         = [VectorField(f.base.space, f.base.name)[f.indices[0]] for f in indexed]
-
-            expr = self.subs(zip(indexed, new_indexed))
-            expr = expr.subs(zip(vector_fields, new_vector_fields))
-            expr = expr.subs(zip(scalar_fields, new_scalar_fields))
-
-        elif self.is_functional:
-            domain = self.domain
-            expr   = self.expr
-            fields = self.fields
-            scalar_fields   = [f for f in fields if isinstance(f.space, ScalarFunctionSpace)]
-            vector_fields   = [f for f in fields if isinstance(f.space, VectorFunctionSpace)]
-
-            new_scalar_fields   = [f.space.field(f.name) for f in scalar_fields]
-            new_vector_fields   = [f.space.field(f.name) for f in vector_fields]
-
-            indexed             = self.expr.atoms(IndexedTestTrial)
-            indexed             = [f for f in indexed if f.base in vector_fields]
-            new_indexed         = [VectorField(f.base.space, f.base.name)[f.indices[0]] for f in indexed]
-
-            expr = expr.subs(zip(indexed, new_indexed))
-            expr = expr.subs(zip(vector_fields, new_vector_fields))
-            expr = expr.subs(zip(scalar_fields, new_scalar_fields))
-            expr = self.func(expr, self.domain, evaluate=False)
-            if self.is_norm:
-                expr._exponent = self._exponent
-
-        expr.is_annotated = True
-        return expr
-
 #==============================================================================
 # TODO check unicity of domain in __new__
 class BasicForm(Expr):
     is_Function = True
     is_linear   = False
     is_bilinear = False
     is_functional = False
-    is_annotated  = False
     is_norm       = False
     _domain     = None
     _ldim        = None
@@ -204,48 +156,3 @@ def _update_free_variables(self, **kwargs):
         # ...
 
         return expr
-
-    def annotate(self):
-
-        if self.is_annotated:
-            return self
-
-        if self.is_bilinear or self.is_linear:
-            fields          = self.fields
-            scalar_fields   = [f for f in fields if isinstance(f.space, ScalarFunctionSpace)]
-            vector_fields   = [f for f in fields if isinstance(f.space, VectorFunctionSpace)]
-
-            new_scalar_fields   = [f.space.field(f.name) for f in scalar_fields]
-            new_vector_fields   = [f.space.field(f.name) for f in vector_fields]
-
-            indexed             = self.expr.atoms(IndexedTestTrial)
-            indexed             = [f for f in indexed if f.base in vector_fields]
-            new_indexed         = [VectorField(f.base.space, f.base.name)[f.indices[0]] for f in indexed]
-
-            expr = self.subs(zip(indexed, new_indexed))
-            expr = expr.subs(zip(vector_fields, new_vector_fields))
-            expr = expr.subs(zip(scalar_fields, new_scalar_fields))
-
-        elif self.is_functional:
-            domain = self.domain
-            expr   = self.expr
-            fields = self.fields
-            scalar_fields   = [f for f in fields if isinstance(f.space, ScalarFunctionSpace)]
-            vector_fields   = [f for f in fields if isinstance(f.space, VectorFunctionSpace)]
-
-            new_scalar_fields   = [f.space.field(f.name) for f in scalar_fields]
-            new_vector_fields   = [f.space.field(f.name) for f in vector_fields]
-
-            indexed             = self.expr.atoms(IndexedTestTrial)
-            indexed             = [f for f in indexed if f.base in vector_fields]
-            new_indexed         = [VectorField(f.base.space, f.base.name)[f.indices[0]] for f in indexed]
-
-            expr = expr.subs(zip(indexed, new_indexed))
-            expr = expr.subs(zip(vector_fields, new_vector_fields))
-            expr = expr.subs(zip(scalar_fields, new_scalar_fields))
-            expr = self.func(expr, self.domain, evaluate=False)
-            if self.is_norm:
-                expr._exponent = self._exponent
-
-        expr.is_annotated = True
-        return expr

sympde/expr/evaluation.py

@@ -399,7 +399,6 @@ def __new__(cls, *args, **options):
         # (Try to) sympify args first
 
         if options.pop('evaluate', True):
-            args = cls._annotate(*args)
             r = cls.eval(*args, **options)
         else:
             r = None
@@ -416,26 +415,6 @@ def __getitem__(self, indices, **kw_args):
         else:
             return Indexed(self, indices, **kw_args)
 
-    # TODO should we keep it?
-    def _annotate(*args):
-        args = list(args)
-        expr = args[0]
-        if isinstance(expr, BasicForm):
-            if not expr.is_annotated:
-                expr = expr.annotate()
-        else:
-            if isinstance(expr, (Add, Mul)):
-                indexed_fields = list(expr.atoms(IndexedTestTrial))
-                new_indexed_fields = [VectorField(F.base.space,F.base.name) for F in indexed_fields]
-                new_indexed_fields = [new_F[F.indices[0]] for new_F,F in zip(new_indexed_fields, indexed_fields)]
-                expr = expr.subs(zip(indexed_fields, new_indexed_fields))
-                fields = list(expr.atoms(ScalarTestFunction,VectorTestFunction).difference(indexed_fields))
-                new_fields = [f.space.field(f.name) for f in fields]
-                expr = expr.subs(zip(fields, new_fields))
-
-        args[0] = expr
-        return args
-
     @classmethod
     @cacheit
     def eval(cls, *_args, **kwargs):
@@ -480,8 +459,10 @@ def eval(cls, *_args, **kwargs):
             J    = expr.mapping.jacobian_expr
             if axis is None:
                 return J
-            else:
+            elif expr.mapping.ldim > 1:
                 return J.col_del(axis)
+            elif expr.mapping.ldim == 1:
+                return J.eye(1)
 
         elif isinstance(expr, JacobianInverseSymbol):
             axis = expr.axis

sympde/topology/basic.py

@@ -199,6 +199,12 @@ def __next__(self):
             raise StopIteration
         self.index += 1
         return result
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        args = ', '.join(sstr(a) for a in self.args)
+        return 'Union({})'.format(args)
+
 #==============================================================================
 class ProductDomain(BasicDomain):
     def __new__(cls, *args, name=None):

sympde/topology/mapping.py

@@ -719,7 +719,6 @@ def eval(cls, F, v):
 #==============================================================================
 class LogicalExpr(CalculusFunction):
 
-    @cacheit
     def __new__(cls, expr, **options):
         # (Try to) sympify args first
 
@@ -756,7 +755,6 @@ def __getitem__(self, indices, **kw_args):
             return Indexed(self, indices, **kw_args)
 
     @classmethod
-    @cacheit
     def eval(cls, expr, mapping=None, dim=None, **options):
         """."""
 
@@ -919,7 +917,7 @@ def eval(cls, expr, mapping=None, dim=None, **options):
 #            v   = cls.eval(grad(expr.args[0]), mapping=mapping, dim=dim)
 #            v   = mapping.jacobian.inv().T*grad(v)
 #            return v
-          
+
         elif isinstance(expr, (dot, inner, outer)):
             args = [cls.eval(arg, mapping=mapping, dim=dim) for arg in expr.args]
             return type(expr)(*args)
@@ -935,9 +933,7 @@ def eval(cls, expr, mapping=None, dim=None, **options):
                 line = []
                 for i_col in range(0, n_cols):
                     line.append(cls.eval(expr[i_row,i_col], mapping=mapping, dim=dim))
-
                 lines.append(line)
-
             return type(expr)(lines)
 
         elif isinstance(expr, dx):
@@ -1028,6 +1024,7 @@ def eval(cls, expr, mapping=None, dim=None, **options):
 
         elif isinstance(expr, (Symbol, Indexed)):
             return expr
+
         elif isinstance(expr, NormalVector):
             return expr
 
@@ -1065,36 +1062,44 @@ def eval(cls, expr, mapping=None, dim=None, **options):
             return Integral(body, domain)
 
         elif isinstance(expr, BilinearForm):
-            tests  = [get_logical_test_function(a) for a in expr.test_functions]
-            trials = [get_logical_test_function(a) for a in expr.trial_functions]
-            body   = cls.eval(expr.expr, mapping=mapping, dim=dim)
+            tests   = [get_logical_test_function(a) for a in expr.test_functions]
+            trials  = [get_logical_test_function(a) for a in expr.trial_functions]
+            mapping = expr.domain.mapping
+            dim     = expr.domain.dim
+            body    = cls.eval(expr.expr, mapping=mapping, dim=dim)
             return BilinearForm((trials, tests), body)
 
         elif isinstance(expr, LinearForm):
-            tests  = [get_logical_test_function(a) for a in expr.test_functions]
-            body   = cls.eval(expr.expr, mapping=mapping, dim=dim)
+            tests   = [get_logical_test_function(a) for a in expr.test_functions]
+            mapping = expr.domain.mapping
+            dim     = expr.domain.dim
+            body    = cls.eval(expr.expr, mapping=mapping, dim=dim)
             return LinearForm(tests, body)
+
         elif isinstance(expr, Norm):
             kind           = expr.kind
             domain         = expr.domain.logical_domain
+            mapping        = expr.domain.mapping
+            dim            = expr.domain.dim
             exponent       = expr.exponent
             e              = cls.eval(expr.expr, mapping=mapping, dim=dim)
             norm           = Norm(e, domain, kind, evaluate=False)
             norm._exponent = exponent
             return norm
+
         elif isinstance(expr, DomainExpression):
             domain  = expr.target.logical_domain
             mapping = expr.target.mapping
             dim     = domain.dim
             J       = mapping.jacobian
             newexpr = cls.eval(expr.expr, mapping=mapping, dim=dim)
             det     = TerminalExpr(sqrt((J.T*J).det()), dim=dim, logical=True, mapping=mapping)
-
             return DomainExpression(domain, ImmutableDenseMatrix([[newexpr*det]]))
 
         elif isinstance(expr, Function):
             args = [cls.eval(a, mapping=mapping, dim=dim) for a in expr.args]
             return type(expr)(*args)
+
         return cls(expr, mapping=mapping, dim=dim, evaluate=False)
 
 #==============================================================================

sympde/topology/space.py

@@ -635,7 +635,7 @@ class ScalarField(Symbol):
     Examples
 
     """
-    _space = None
+    _space         = None
     is_commutative = True
     _projection_of = None
     def __new__(cls, space, name=None):
@@ -661,6 +661,14 @@ def projection_of(self):
     def set_as_projection(self, expr):
         self._projection_of = expr
 
+    def __eq__(self, a):
+        if isinstance(a, ScalarField):
+            return a.space == self.space and a.name == self.name
+        return False
+
+    def __hash__(self):
+        return hash((self.name, self.space))
+
     def _sympystr(self, printer):
         sstr = printer.doprint
         return sstr(self.name)

sympde/version.py

@@ -1 +1 @@
-__version__ = "0.10.4"
+__version__ = "0.10.5"