From e1df7a9f8455ae7e1b7bc162a3c27bd81c5cf415 Mon Sep 17 00:00:00 2001 From: Said Hadjout Date: Thu, 21 Jan 2021 16:36:21 +0000 Subject: [PATCH] Feec dirichlet (#87) The following changes were made: - Delete node annotation; - Compute boundary jacobian when dim = 1; - Update library version. --- sympde/expr/basic.py | 93 -------------------------------------- sympde/expr/evaluation.py | 25 ++-------- sympde/topology/basic.py | 6 +++ sympde/topology/mapping.py | 27 ++++++----- sympde/topology/space.py | 10 +++- sympde/version.py | 2 +- 6 files changed, 35 insertions(+), 128 deletions(-) diff --git a/sympde/expr/basic.py b/sympde/expr/basic.py index 465105c5..4cc9e407 100644 --- a/sympde/expr/basic.py +++ b/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,52 +84,6 @@ 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): @@ -138,7 +91,6 @@ class BasicForm(Expr): 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 diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index fa955e04..95ba51f3 100644 --- a/sympde/expr/evaluation.py +++ b/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 diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py index a1dae1de..517aebca 100644 --- a/sympde/topology/basic.py +++ b/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): diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index f87729ba..ed9778c5 100644 --- a/sympde/topology/mapping.py +++ b/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,23 +1062,31 @@ 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 @@ -1089,12 +1094,12 @@ def eval(cls, expr, mapping=None, dim=None, **options): 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) #============================================================================== diff --git a/sympde/topology/space.py b/sympde/topology/space.py index 3a5ddf82..bb89634a 100644 --- a/sympde/topology/space.py +++ b/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) diff --git a/sympde/version.py b/sympde/version.py index d9b054ab..a67aac09 100644 --- a/sympde/version.py +++ b/sympde/version.py @@ -1 +1 @@ -__version__ = "0.10.4" +__version__ = "0.10.5"