From 9445c96acaf3d65b89f3723163e99df042d1961e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 18 Apr 2023 11:20:48 +0200 Subject: [PATCH 1/4] Add PushForward class to sympde.topology.mapping --- sympde/topology/mapping.py | 150 +++++++++++++++++++++++++++++++------ 1 file changed, 129 insertions(+), 21 deletions(-) diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 5c8b6589..4528394d 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -858,6 +858,74 @@ def eval(cls, F, v): v = M*v return Tuple(*v) +#============================================================================== +def get_kind_type(kind): + + from .datatype import SpaceType, dtype_space_registry + + if kind is None: + kind = 'undefined' + + if isinstance(kind, str): + kind_str = kind.lower() + assert kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined'] + kind_type = dtype_space_registry[kind_str] + + elif isinstance(kind, SpaceType): + kind_type = kind + + else: + raise TypeError('Expecting kind to be string or SpaceType') + + return kind_type + +#============================================================================== +class PushForward(Expr): + is_commutative = False + + def __new__(cls, expr, *, mapping=None, kind='h1', evaluate=False): + + kind = get_kind_type(kind) + + if evaluate: + return cls.eval(expr, mapping=mapping, kind=kind) + + obj = Expr.__new__(cls, expr, mapping, kind) + return obj + + @property + def arg(self): + return self.args[0] + + @property + def mapping(self): + return self.args[1] + + @property + def kind(self): + return self.args[2] + + @classmethod + def eval(cls, expr, mapping, kind): + + kind = get_kind_type(kind) + + if kind in ('h1', None): + return expr + + J = mapping.jacobian + + if kind is HcurlSpaceType(): + return J.inv().T * expr + + if kind is HdivSpaceType(): + return (J / J.det()) * expr + + if kind is L2SpaceType(): + return expr / J.det() + + raise ValueError(f'Cannot compute Push-Forward of kind = {kind}') + #============================================================================== class LogicalExpr(CalculusFunction): @@ -898,7 +966,9 @@ def eval(cls, expr, domain, **options): from sympde.expr.expr import BilinearForm, LinearForm, BasicForm, Norm from sympde.expr.expr import Integral - types = (ScalarFunction, VectorFunction, DifferentialOperator, Trace, Integral) + # Expression types for which special rules are applied + types = (ScalarFunction, VectorFunction, DifferentialOperator, + Trace, Integral, PushForward) mapping = domain.mapping dim = domain.dim @@ -971,12 +1041,25 @@ def eval(cls, expr, domain, **options): return newexpr elif isinstance(expr, (VectorFunction, ScalarFunction)): - return PullBack(expr, mapping).expr + u = expr + u_hat = get_logical_test_function(u) + space = u.space + mapping = space.domain.mapping + kind = space.kind + return PushForward(u_hat, mapping=mapping, kind=kind, evaluate=True) + +# elif isinstance(expr, (VectorFunction, ScalarFunction)): +# return PullBack(expr, mapping).expr elif isinstance(expr, Transpose): arg = cls(expr.arg, domain) return Transpose(arg) - + + elif isinstance(expr, PushForward): + # Remove push-forward operation + arg = expr.args[0] + return cls(arg, domain) + elif isinstance(expr, grad): arg = expr.args[0] if isinstance(mapping, InterfaceMapping): @@ -989,40 +1072,65 @@ def eval(cls, expr, domain, **options): else: raise TypeError(arg) - arg = type(arg)(cls.eval(a, domain)) + arg_log = type(arg)(cls.eval(a, domain)) else: - arg = cls.eval(arg, domain) + arg_log = cls.eval(arg, domain) - return mapping.jacobian.inv().T*grad(arg) + return PushForward(grad(arg_log), mapping=mapping, kind='hcurl', evaluate=True) +# return mapping.jacobian.inv().T * grad(arg_log) elif isinstance(expr, curl): + arg = expr.args[0] + + # If on interface, select which mapping to use if isinstance(mapping, InterfaceMapping): if isinstance(arg, MinusInterfaceOperator): - arg = arg.args[0] + int_arg = arg.args[0] mapping = mapping.minus elif isinstance(arg, PlusInterfaceOperator): - arg = arg.args[0] + int_arg = arg.args[0] mapping = mapping.plus else: raise TypeError(arg) - if isinstance(arg, VectorFunction): - arg = PullBack(arg, mapping) + J = mapping.jacobian + int_arg_log = cls.eval(int_arg, domain) + arg_log = type(expr.args[0])(int_arg_log) + else: - arg = cls.eval(arg, domain) + arg_log = cls.eval(arg, domain) - if isinstance(arg, PullBack) and isinstance(arg.kind, HcurlSpaceType): - J = mapping.jacobian - arg = arg.test - if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)): - arg = type(expr.args[0])(arg) - if expr.is_scalar: - return (1/J.det())*curl(arg) + return PushForward(curl(arg_log), mapping=mapping, kind='hdiv', evaluate=True) - return (J/J.det())*curl(arg) - else: - raise NotImplementedError('TODO') +# elif isinstance(expr, curl): +# arg = expr.args[0] +# if isinstance(mapping, InterfaceMapping): +# if isinstance(arg, MinusInterfaceOperator): +# arg = arg.args[0] +# mapping = mapping.minus +# elif isinstance(arg, PlusInterfaceOperator): +# arg = arg.args[0] +# mapping = mapping.plus +# else: +# raise TypeError(arg) +# +# if isinstance(arg, VectorFunction): +# arg = PullBack(arg, mapping) +# else: +# arg = cls.eval(arg, domain) +# +# if isinstance(arg, PullBack) and isinstance(arg.kind, HcurlSpaceType): +# J = mapping.jacobian +# arg = arg.test +# if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)): +# arg = type(expr.args[0])(arg) +# if expr.is_scalar: +# return (1/J.det())*curl(arg) +# +# return (J/J.det())*curl(arg) +# else: +# raise NotImplementedError('TODO') elif isinstance(expr, div): arg = expr.args[0] From 9a82f4e74abb82cc16dff93385b4b99e91c3f29b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 18 Apr 2023 11:46:51 +0200 Subject: [PATCH 2/4] Compare kind to H1SpaceType() instead of 'h1' string --- sympde/topology/mapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 4528394d..fc83e54b 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -910,7 +910,7 @@ def eval(cls, expr, mapping, kind): kind = get_kind_type(kind) - if kind in ('h1', None): + if kind in (H1SpaceType(), None): return expr J = mapping.jacobian From 2fb6306b2c6505cbf4c1b53cfa1b12b3cbfb3df1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 18 Apr 2023 11:52:35 +0200 Subject: [PATCH 3/4] Compare kind to UndefinedSpaceType() instead of None --- sympde/topology/mapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index fc83e54b..34fdb9b5 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -910,7 +910,7 @@ def eval(cls, expr, mapping, kind): kind = get_kind_type(kind) - if kind in (H1SpaceType(), None): + if kind in (H1SpaceType(), UndefinedSpaceType()): return expr J = mapping.jacobian From f7f8b82b4d6d0b4718476e48f26d04205718bdbf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 21 Apr 2023 13:50:05 +0200 Subject: [PATCH 4/4] Add new class PullBack, which calls LogicalExpr if needed --- sympde/topology/mapping.py | 239 ++++++++++++++++++++++++------------- 1 file changed, 159 insertions(+), 80 deletions(-) diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 34fdb9b5..339df5b0 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -879,6 +879,89 @@ def get_kind_type(kind): return kind_type +#============================================================================== +class PullBack(Expr): + is_commutative = False + + def __new__(cls, expr, *, mapping=None, domain=None, kind='h1', evaluate=False): + + kind = get_kind_type(kind) + + if evaluate: + return cls.eval(expr, mapping=mapping, domain=domain, kind=kind) + + obj = Expr.__new__(cls, expr, mapping, kind) + return obj + + @property + def expr(self): + return self.args[0] + + @property + def mapping(self): + return self.args[1] + + @property + def kind(self): + return self.args[2] + + @classmethod + def eval(cls, expr, domain, mapping, kind): + + kind = get_kind_type(kind) + + if isinstance(expr, (ScalarFunction, VectorFunction)): + u = expr + u_hat = get_logical_test_function(u) + space = u.space + mapping = space.domain.mapping + + # By definition u is the PushForward of u_hat. (This is implicit.) + # Hence, if the PullBack and the PushForward are of the same kind, + # they cancel each other. + if kind == u.space.kind: + return u_hat + + # If the kinds do not match, continue processing... + + if isinstance(expr, PushForward): + # If PullBack and PushForward are of the same kind, they simplify. + if kind == expr.kind: + return expr.arg + + # Otherwise, evaluate PushForward explicitly and redefine expression + else: + expr = PushForward( + expr.arg, + mapping = expr.mapping, + kind = expr.kind, + evaluate = True + ) + + # In all other cases use LogicalExpr.eval() and apply inverse formulas + # TODO: get domain + log_expr = LogicalExpr.eval(expr, domain) + + # TODO: verify what this means + if space.is_broken: + assert mapping is not None + + if kind in (H1SpaceType(), UndefinedSpaceType()): + return log_expr + + J = mapping.jacobian + + if kind is HcurlSpaceType(): + return J.T * log_expr + + if kind is HdivSpaceType(): + return J.det() * J.inv() * log_expr + + if kind is L2SpaceType(): + return J.det() * log_expr + + raise ValueError(f'Cannot compute Pull-Back of kind = {kind}') + #============================================================================== class PushForward(Expr): is_commutative = False @@ -980,6 +1063,7 @@ def eval(cls, expr, domain, **options): if not has(expr, types): if has(expr, DiffOperator): + # TODO: do not rewrite return cls( expr, domain, evaluate=False) else: syms = symbols(ph_coords[:dim]) @@ -1040,69 +1124,64 @@ def eval(cls, expr, domain, **options): newexpr = newexpr.expr.subs(test, PlusInterfaceOperator(test)) return newexpr +# elif isinstance(expr, (VectorFunction, ScalarFunction)): +# return PullBack(expr, mapping).expr + elif isinstance(expr, (VectorFunction, ScalarFunction)): u = expr u_hat = get_logical_test_function(u) space = u.space - mapping = space.domain.mapping kind = space.kind + + # TODO: verify what this means + if space.is_broken: + assert mapping is not None + else: + mapping = space.domain.mapping + return PushForward(u_hat, mapping=mapping, kind=kind, evaluate=True) -# elif isinstance(expr, (VectorFunction, ScalarFunction)): -# return PullBack(expr, mapping).expr + elif isinstance(expr, PushForward): + return PushForward(expr.arg, mapping=expr.mapping, kind=expr.kind, + evaluate=True) elif isinstance(expr, Transpose): arg = cls(expr.arg, domain) return Transpose(arg) - elif isinstance(expr, PushForward): - # Remove push-forward operation - arg = expr.args[0] - return cls(arg, domain) - elif isinstance(expr, grad): - arg = expr.args[0] - if isinstance(mapping, InterfaceMapping): - if isinstance(arg, MinusInterfaceOperator): - a = arg.args[0] - mapping = mapping.minus - elif isinstance(arg, PlusInterfaceOperator): - a = arg.args[0] - mapping = mapping.plus - else: - raise TypeError(arg) - - arg_log = type(arg)(cls.eval(a, domain)) - else: - arg_log = cls.eval(arg, domain) - + arg_log = PullBack(expr.args[0], mapping=mapping, domain=domain, kind='h1', evaluate=True) return PushForward(grad(arg_log), mapping=mapping, kind='hcurl', evaluate=True) -# return mapping.jacobian.inv().T * grad(arg_log) elif isinstance(expr, curl): - - arg = expr.args[0] - - # If on interface, select which mapping to use - if isinstance(mapping, InterfaceMapping): - if isinstance(arg, MinusInterfaceOperator): - int_arg = arg.args[0] - mapping = mapping.minus - elif isinstance(arg, PlusInterfaceOperator): - int_arg = arg.args[0] - mapping = mapping.plus - else: - raise TypeError(arg) - - J = mapping.jacobian - int_arg_log = cls.eval(int_arg, domain) - arg_log = type(expr.args[0])(int_arg_log) - - else: - arg_log = cls.eval(arg, domain) - + caller = lambda e, cls=cls: cls.eval(e, cls.domain) + arg_log = PullBack(expr.args[0], mapping=mapping, domain=domain, kind='hcurl', evaluate=True) return PushForward(curl(arg_log), mapping=mapping, kind='hdiv', evaluate=True) + elif isinstance(expr, div): + caller = lambda e, cls=cls: cls.eval(e, cls.domain) + arg_log = PullBack(expr.args[0], mapping=mapping, domain=domain, kind='hdiv', evaluate=True) + return PushForward(div(arg_log), mapping=mapping, kind='l2', evaluate=True) + +# elif isinstance(expr, grad): +# arg = expr.args[0] +# print(f'arg = {arg}') +# if isinstance(mapping, InterfaceMapping): +# if isinstance(arg, MinusInterfaceOperator): +# a = arg.args[0] +# mapping = mapping.minus +# elif isinstance(arg, PlusInterfaceOperator): +# a = arg.args[0] +# mapping = mapping.plus +# else: +# raise TypeError(arg) +# +# arg_log = type(arg)(cls.eval(a, domain)) +# else: +# arg_log = cls.eval(arg, domain) +# +# return mapping.jacobian.inv().T * grad(arg_log) +# # elif isinstance(expr, curl): # arg = expr.args[0] # if isinstance(mapping, InterfaceMapping): @@ -1131,41 +1210,41 @@ def eval(cls, expr, domain, **options): # return (J/J.det())*curl(arg) # else: # raise NotImplementedError('TODO') - - elif isinstance(expr, div): - arg = expr.args[0] - if isinstance(mapping, InterfaceMapping): - if isinstance(arg, MinusInterfaceOperator): - arg = arg.args[0] - mapping = mapping.minus - elif isinstance(arg, PlusInterfaceOperator): - arg = arg.args[0] - mapping = mapping.plus - else: - raise TypeError(arg) - - if isinstance(arg, (ScalarFunction, VectorFunction)): - arg = PullBack(arg, mapping) - else: - - arg = cls.eval(arg, domain) - - if isinstance(arg, PullBack) and isinstance(arg.kind, HdivSpaceType): - J = mapping.jacobian - arg = arg.test - if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)): - arg = type(expr.args[0])(arg) - return (1/J.det())*div(arg) - elif isinstance(arg, PullBack): - return SymbolicTrace(mapping.jacobian.inv().T*grad(arg.test)) - else: - raise NotImplementedError('TODO') - - elif isinstance(expr, laplace): - arg = expr.args[0] - v = cls.eval(grad(arg), domain) - v = mapping.jacobian.inv().T*grad(v) - return SymbolicTrace(v) +# +# elif isinstance(expr, div): +# arg = expr.args[0] +# if isinstance(mapping, InterfaceMapping): +# if isinstance(arg, MinusInterfaceOperator): +# arg = arg.args[0] +# mapping = mapping.minus +# elif isinstance(arg, PlusInterfaceOperator): +# arg = arg.args[0] +# mapping = mapping.plus +# else: +# raise TypeError(arg) +# +# if isinstance(arg, (ScalarFunction, VectorFunction)): +# arg = PullBack(arg, mapping) +# else: +# +# arg = cls.eval(arg, domain) +# +# if isinstance(arg, PullBack) and isinstance(arg.kind, HdivSpaceType): +# J = mapping.jacobian +# arg = arg.test +# if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)): +# arg = type(expr.args[0])(arg) +# return (1/J.det())*div(arg) +# elif isinstance(arg, PullBack): +# return SymbolicTrace(mapping.jacobian.inv().T*grad(arg.test)) +# else: +# raise NotImplementedError('TODO') +# +# elif isinstance(expr, laplace): +# arg = expr.args[0] +# v = cls.eval(grad(arg), domain) +# v = mapping.jacobian.inv().T*grad(v) +# return SymbolicTrace(v) # elif isinstance(expr, hessian): # arg = expr.args[0]