From 29b7ca845c2c2edcad2cd6338e59f7c807fb0ac6 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Fri, 16 Jun 2023 14:53:26 +0200 Subject: [PATCH 1/9] Add codomain_complex flag into ScalarFunctionSpace, VectorFunctionSpace, Derham and create an SesquilinearForm --- sympde/expr/expr.py | 188 +++++++++++++++++++++++++++++++++++++++ sympde/topology/space.py | 36 ++++---- 2 files changed, 208 insertions(+), 16 deletions(-) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 7428df7f..1114ef7e 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -5,6 +5,7 @@ from sympy import Dummy from sympy import Matrix, ImmutableDenseMatrix +from sympy import conjugate from sympy.core import Basic, S from sympy.core import Expr, Add, Mul from sympy.core.numbers import Zero as sy_Zero @@ -430,6 +431,12 @@ def __new__(cls, arguments, expr, **options): # Distinguish between trial and test functions trial_functions, test_functions = args + if trial_functions.space.codomain_complex != test_functions.space.codomain_complex: + raise TypeError('Trial space and Test space should both be real.') + if trial_functions.space.codomain_complex==False: + raise TypeError('Trial space and Test space should be real. In the complex case, a SesquilinearForm should be called ') + + # Check linearity with respect to trial functions if not is_linear_expression(expr, trial_functions): msg = ' Expression is not linear w.r.t trial functions {}'\ @@ -515,6 +522,118 @@ def __call__(self, trials, tests, **kwargs): return expr +#============================================================================== +class SesquilinearForm(BasicForm): + is_sesquilinear = True + _is_hermitian = None + + def __new__(cls, arguments, expr, **options): + + # Trivial case: null expression + if expr == 0: + return sy_Zero() + + # Check that integral expression is given + + if not expr.atoms(Integral): + raise ValueError('Expecting integral Expression') + + # TODO: why do we 'sanitize' here? + args = _sanitize_arguments(arguments, is_bilinear=True) + + # Distinguish between trial and test functions + trial_functions, test_functions = args + if trial_functions.space.codomain_complex != test_functions.space.codomain_complex: + raise TypeError('Trial space and Test space should both be complex.') + if trial_functions.space.codomain_complex==True: + raise TypeError('Trial space and Test space should be complex. In the real case, a BilinearForm should be called ') + + + # Check linearity with respect to trial functions + if not is_antilinear_expression(expr, trial_functions): + msg = ' Expression is not linear w.r.t trial functions {}'\ + .format(trial_functions) + raise UnconsistentLinearExpressionError(msg) + + # Check linearity with respect to test functions + if not is_linear_expression(expr, test_functions): + msg = ' Expression is not linear w.r.t test functions {}'\ + .format(test_functions) + raise UnconsistentLinearExpressionError(msg) + + # Create new object of type SesquilinearForm + obj = Basic.__new__(cls, args, expr) + + # Compute 'domain' property (scalar or tuple) + # TODO: is this is useful? + obj._domain = _get_domain(expr) + + return obj + + @property + def variables(self): + return self._args[0] + + @property + def expr(self): + return self._args[1] + + @property + def domain(self): + return self._domain + + @property + def test_functions(self): + return self.variables[1] + + @property + def trial_functions(self): + return self.variables[0] + + @property + def test_spaces(self): + return [u.space for u in self.test_functions] + + @property + def trial_spaces(self): + return [u.space for u in self.trial_functions] + + @property + def coordinates(self): + return self.test_spaces[0].coordinates + + @property + def ldim(self): + return self.test_spaces[0].ldim + + @property + def is_hermitian(self): + if self._is_hermitian is None: + left, right = self.variables + a1 = self(left, right) + a2 = self(right, left) + self._is_hermitian = (a1 == a2) + return self._is_hermitian + + def __call__(self, trials, tests, **kwargs): + + # Use free variables if given and available + expr = self._update_free_variables(**kwargs) + + # If needed, convert positional arguments to lists + if not is_sequence(trials): trials = [trials] + if not is_sequence(tests ): tests = [tests ] + + # Concatenate input values into single list + values = [*trials, *tests] + + # Substitute variables with given values in bilinear expression + variables = [*self.variables[0], *self.variables[1]] + subs = dict(zip(variables, values)) + expr, _ = expr._xreplace(subs) + + return expr + #============================================================================== class Norm(Functional): is_norm = True @@ -733,6 +852,75 @@ def is_linear_expression(expr, args, integral=True, debug=True): return True +#============================================================================== +def is_antilinear_expression(expr, args, debug=True): + """checks if an expression is linear with respect to the given arguments.""" + # ... + left_args = [] + right_args = [] + + for arg in args: + tag = random_string( 4 ) + + if isinstance(arg, ScalarFunction): + left = ScalarFunction(arg.space, name='l_' + tag) + right = ScalarFunction(arg.space, name='r_' + tag) + + elif isinstance(arg, VectorFunction): + left = VectorFunction(arg.space, name='l_' + tag) + right = VectorFunction(arg.space, name='r_' + tag) + else: + raise TypeError('argument must be a {Scalar|Vector}Function') + + left_args += [left] + right_args += [right] + # ... + + # ... check addition + newargs = [left + right for left, right in zip(left_args, right_args)] + + newexpr = expr.subs(zip(args, newargs)) + left_expr = expr.subs(zip(args, left_args)) + right_expr = expr.subs(zip(args, right_args)) + + a = newexpr + b = left_expr + right_expr + + if not( (a-b).expand() == 0 or a.expand() == b.expand()): + # TODO use a warning or exception? + if debug: + print('Failed to assert addition property') + print('{} != {}'.format(a.expand(), b.expand())) + return False + + # ... + + # ... check multiplication + tag = random_string( 4 ) + coeff = Constant('alpha_' + tag) + + newexpr = expr + for arg, left in zip(args, left_args): + newarg = coeff * left + newexpr = newexpr.subs(arg, newarg) + + atoms = list(newexpr.atoms(BasicOperator)) + subs = [e.func(*e.args, evaluate=True) for e in atoms] + newexpr = newexpr.subs(zip(atoms, subs)) + + + left_expr = expr.subs(list(zip(args, left_args))) + left_expr = conjugate(coeff) * left_expr + if not( (newexpr-left_expr).expand() == 0 or newexpr.expand()==left_expr.expand()): + # TODO use a warning or exception? + if debug: + print('Failed to assert multiplication property') + print('{} != {}'.format(newexpr, left_expr)) + return False + # ... + + return True + #============================================================================== def integral(domain, expr): return Integral(expr, domain) diff --git a/sympde/topology/space.py b/sympde/topology/space.py index cce1f556..96dffa7d 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -135,7 +135,7 @@ class BasicFunctionSpace(Basic): _kind = None _regularity = None # TODO pass it as arg to __new__ _is_broken = None - def __new__(cls, name, domain, shape, kind): + def __new__(cls, name, domain, shape, kind, *, codomain_complex=False): if not isinstance(domain, BasicDomain): raise TypeError('> Expecting a BasicDomain object for domain') @@ -144,7 +144,7 @@ def __new__(cls, name, domain, shape, kind): obj._name = name obj._domain = domain obj._shape = shape - + obj._codomain_complex=codomain_complex # ... if kind is None: kind = 'undefined' @@ -206,6 +206,10 @@ def is_broken(self): def regularity(self): return self._regularity + @property + def codomain_type(self): + return self._codomain_type + @property def coordinates(self): return self.domain.coordinates @@ -224,9 +228,9 @@ class ScalarFunctionSpace(BasicFunctionSpace): """ Represents a basic continuous scalar Function space. """ - def __new__(cls, name, domain, kind=None): + def __new__(cls, name, domain, kind=None, *, codomain_complex=False): shape = 1 - return BasicFunctionSpace.__new__(cls, name, domain, shape, kind) + return BasicFunctionSpace.__new__(cls, name, domain, shape, kind, codomain_complex=codomain_complex) def element(self, name): return ScalarFunction(self, name) @@ -236,9 +240,9 @@ class VectorFunctionSpace(BasicFunctionSpace): """ Represents a basic continuous vector Function space. """ - def __new__(cls, name, domain, kind=None): + def __new__(cls, name, domain, kind=None, *, codomain_complex=False): shape = domain.dim - return BasicFunctionSpace.__new__(cls, name, domain, shape, kind) + return BasicFunctionSpace.__new__(cls, name, domain, shape, kind, codomain_complex=codomain_complex) def element(self, name): return VectorFunction(self, name) @@ -246,7 +250,7 @@ def element(self, name): #============================================================================= class Derham: """.""" - def __init__(self, domain, sequence=None): + def __init__(self, domain, sequence=None, *, codomain_complex=False): shape = domain.dim self._V0 = None self._V1 = None @@ -254,8 +258,8 @@ def __init__(self, domain, sequence=None): self._V3 = None if shape == 1: - spaces = [ScalarFunctionSpace('H1', domain, kind='H1'), - ScalarFunctionSpace('L2', domain, kind='L2')] + spaces = [ScalarFunctionSpace('H1', domain, kind='H1', codomain_complex=codomain_complex), + ScalarFunctionSpace('L2', domain, kind='L2', codomain_complex=codomain_complex)] self._V0 = spaces[0] self._V1 = spaces[1] @@ -264,19 +268,19 @@ def __init__(self, domain, sequence=None): assert sequence is not None space = sequence[1] - spaces = [ScalarFunctionSpace('H1', domain, kind='H1'), - VectorFunctionSpace(space, domain, kind=space), - ScalarFunctionSpace('L2', domain, kind='L2')] + spaces = [ScalarFunctionSpace('H1', domain, kind='H1' ,codomain_complex=codomain_complex), + VectorFunctionSpace(space, domain, kind=space, codomain_complex=codomain_complex), + ScalarFunctionSpace('L2', domain, kind='L2', codomain_complex=codomain_complex)] self._V0 = spaces[0] self._V1 = spaces[1] self._V2 = spaces[2] elif shape == 3: - spaces = [ScalarFunctionSpace('H1', domain, kind='H1'), - VectorFunctionSpace('Hcurl', domain, kind='Hcurl'), - VectorFunctionSpace('Hdiv', domain, kind='Hdiv'), - ScalarFunctionSpace('L2', domain, kind='L2')] + spaces = [ScalarFunctionSpace('H1', domain, kind='H1' , codomain_complex=codomain_complex), + VectorFunctionSpace('Hcurl', domain, kind='Hcurl', codomain_complex=codomain_complex), + VectorFunctionSpace('Hdiv', domain, kind='Hdiv' , codomain_complex=codomain_complex), + ScalarFunctionSpace('L2', domain, kind='L2' , codomain_complex=codomain_complex)] self._V0 = spaces[0] self._V1 = spaces[1] From a99f20c6012055425cd0733420f8fac3fad3d980 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Mon, 19 Jun 2023 14:43:07 +0200 Subject: [PATCH 2/9] Make the dot product sustain 1D product, correct the test in BilinearForm --- sympde/calculus/core.py | 13 ++++++++++--- sympde/expr/expr.py | 4 ++-- sympde/topology/space.py | 10 +++++----- 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index ac831aaa..2bb25312 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -96,7 +96,7 @@ from sympy import Indexed, sympify from sympy import Matrix, ImmutableDenseMatrix -from sympy import cacheit +from sympy import cacheit, conjugate from sympy.core import Basic from sympy.core import Add, Mul, Pow from sympy.core.containers import Tuple @@ -104,6 +104,8 @@ from sympy.core.decorators import call_highest_priority from sympy.core.compatibility import is_sequence +from math import sqrt + from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery from sympde.topology.space import ScalarFunction, VectorFunction @@ -266,6 +268,7 @@ def is_scalar(atom): """ return is_constant(atom) or isinstance(atom, ScalarFunction) + #============================================================================== # TODO add dot(u,u) +2*dot(u,v) + dot(v,v) = dot(u+v,u+v) # now we only have dot(u,u) + dot(u,v)+ dot(v,u) + dot(v,v) = dot(u+v,u+v) @@ -351,9 +354,14 @@ def __new__(cls, arg1, arg2, **options): args_2 = [i for i in b if not i.is_commutative] c2 = [i for i in b if not i in args_2] + + c = Mul(*c1)*Mul(*c2) + # 1D case where everything is commutative + if args_1==[] and args_2==[]: + return(c) + a = reduce(mul, args_1) b = reduce(mul, args_2) - c = Mul(*c1)*Mul(*c2) if str(a) > str(b): a,b = b,a @@ -780,7 +788,6 @@ def eval(cls, expr): raise ArgumentTypeError(msg) return cls(expr, evaluate=False) - #============================================================================== class Curl(DiffOperator): """ diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 1114ef7e..615eb813 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -431,9 +431,9 @@ def __new__(cls, arguments, expr, **options): # Distinguish between trial and test functions trial_functions, test_functions = args - if trial_functions.space.codomain_complex != test_functions.space.codomain_complex: + if trial_functions[0].space.codomain_complex != test_functions[0].space.codomain_complex: raise TypeError('Trial space and Test space should both be real.') - if trial_functions.space.codomain_complex==False: + if trial_functions[0].space.codomain_complex==True: raise TypeError('Trial space and Test space should be real. In the complex case, a SesquilinearForm should be called ') diff --git a/sympde/topology/space.py b/sympde/topology/space.py index 96dffa7d..3f829b72 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -207,8 +207,8 @@ def regularity(self): return self._regularity @property - def codomain_type(self): - return self._codomain_type + def codomain_complex(self): + return self._codomain_complex @property def coordinates(self): @@ -442,9 +442,9 @@ class ScalarFunction(Symbol): def __new__(cls, space, name): if not isinstance(space, ScalarFunctionSpace): raise ValueError('Expecting a ScalarFunctionSpace') - obj = Expr.__new__(cls) - obj._space = space - obj._name = name + obj = Expr.__new__(cls) + obj._space = space + obj._name = name return obj @property From 6fc438ee813596cb563e9403b99650b768a905e7 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Tue, 20 Jun 2023 09:13:25 +0200 Subject: [PATCH 3/9] Allow the class TerminalExpr and DifferentialOperator to read the expression conjugate(), Add the SesquilinearForm to the importable function, --- sympde/expr/evaluation.py | 6 +++++- sympde/expr/expr.py | 1 + sympde/topology/derivatives.py | 5 ++++- 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 7cd02bbe..3b1fa1c8 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -6,7 +6,7 @@ import numpy as np from sympy import Abs, S, cacheit -from sympy import Indexed, Matrix, ImmutableDenseMatrix +from sympy import Indexed, Matrix, ImmutableDenseMatrix, conjugate from sympy import expand from sympy.core import Basic, Symbol from sympy.core import Add, Mul, Pow @@ -828,6 +828,10 @@ def eval(cls, expr, domain): expr = cls(expr.expr, domain=domain) return LogicalExpr(expr, domain=domain) + elif isinstance(expr, conjugate): + expr = cls(expr.args[0], domain=domain) + return conjugate(expr) + return expr diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 615eb813..34e29f48 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -31,6 +31,7 @@ __all__ = ( 'BilinearForm', + 'SesquilinearForm', 'Functional', 'IntAdd', 'Integral', diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index b1b6ed0a..f6cdc688 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -15,7 +15,7 @@ from sympy import S from sympy import Indexed from sympy import diff -from sympy import log +from sympy import log, conjugate from sympy import preorder_traversal from sympy import cacheit from sympy.core.compatibility import is_sequence @@ -162,6 +162,9 @@ def eval(cls, expr): v = (log(b)*cls(e, evaluate=True) + e*cls(b, evaluate=True)/b) * b**e return v + elif isinstance(expr, conjugate): + return conjugate(cls(expr.args[0], evaluate=True)) + else: msg = '{expr} of type {type}'.format(expr=expr, type=type(expr)) raise NotImplementedError(msg) From e657d25505c5971160999412eabfb5f2880531d5 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Tue, 20 Jun 2023 10:53:20 +0200 Subject: [PATCH 4/9] Adapt BasicExpr, Equation, LogicalExpr for SesquilinearForm. Correct assertion on SesquilinearForm --- sympde/expr/basic.py | 6 ++++-- sympde/expr/equation.py | 11 +++++++---- sympde/expr/evaluation.py | 2 +- sympde/expr/expr.py | 11 +++++++---- sympde/topology/mapping.py | 10 ++++++++-- 5 files changed, 27 insertions(+), 13 deletions(-) diff --git a/sympde/expr/basic.py b/sympde/expr/basic.py index 96660966..4db54cfb 100644 --- a/sympde/expr/basic.py +++ b/sympde/expr/basic.py @@ -62,14 +62,16 @@ class BasicExpr(Expr): is_Function = True is_linear = False is_bilinear = False + # TODO put is_sesquilinear=False here but it don't work + #is_sesquilinear = False is_functional = False @property def fields(self): atoms = self.expr.atoms(ScalarFunction, VectorFunction) - if self.is_bilinear or self.is_linear: + if self.is_bilinear or self.is_sesquilinear or self.is_linear: args = self.variables - if self.is_bilinear: + if self.is_bilinear or self.is_sesquilinear : args = args[0]+args[1] fields = tuple(atoms.difference(args)) else: diff --git a/sympde/expr/equation.py b/sympde/expr/equation.py index 5b4ed1aa..137043d5 100644 --- a/sympde/expr/equation.py +++ b/sympde/expr/equation.py @@ -14,7 +14,7 @@ from sympde.calculus import grad, dot from sympde.core.utils import random_string -from .expr import BilinearForm, LinearForm +from .expr import BilinearForm, LinearForm, SesquilinearForm from .expr import linearize from .errors import ( UnconsistentLhsError, UnconsistentRhsError, UnconsistentArgumentsError, UnconsistentBCError ) @@ -192,8 +192,8 @@ class Equation(Basic): def __new__(cls, lhs, rhs, trials, tests, bc=None, constraint=None): # ... - if not isinstance(lhs, BilinearForm): - raise UnconsistentLhsError('> lhs must be a bilinear') + if not isinstance(lhs, (BilinearForm, SesquilinearForm)): + raise UnconsistentLhsError('> lhs must be a bilinear or a Sesquilinear') if not isinstance(rhs, LinearForm): raise UnconsistentRhsError('> rhs must be a linear') @@ -380,7 +380,10 @@ def __new__(cls, form, fields, bc=None, trials=None): def find(trials, *, forall, lhs, rhs, bc=None, constraint=None): tests = forall - lhs = BilinearForm((trials, tests), lhs) + if tests.space.codomain_complex: + lhs = SesquilinearForm((trials, tests), lhs) + else: + lhs = BilinearForm((trials, tests), lhs) rhs = LinearForm( tests , rhs) return Equation(lhs, rhs, trials, tests, bc=bc, constraint=constraint) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 3b1fa1c8..6516b28d 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -139,7 +139,7 @@ def _get_trials_tests(expr, *, flatten=False): if not isinstance(expr, (BasicForm, BasicExpr)): raise TypeError("Expression must be of type BasicForm or BasicExpr, got '{}' instead".format(type(expr))) - if expr.is_bilinear: + if expr.is_bilinear or expr.is_sesquilinear: trials = _unpack_functions(expr.variables[0]) if flatten else expr.variables[0] tests = _unpack_functions(expr.variables[1]) if flatten else expr.variables[1] diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 34e29f48..b9a949c2 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -334,6 +334,7 @@ def space(self): #============================================================================== class LinearForm(BasicForm): is_linear = True + is_sesquilinear = False def __new__(cls, arguments, expr, **options): @@ -413,6 +414,7 @@ def __call__(self, *tests, **kwargs): #============================================================================== class BilinearForm(BasicForm): is_bilinear = True + is_sesquilinear = False _is_symmetric = None def __new__(cls, arguments, expr, **options): @@ -434,7 +436,7 @@ def __new__(cls, arguments, expr, **options): if trial_functions[0].space.codomain_complex != test_functions[0].space.codomain_complex: raise TypeError('Trial space and Test space should both be real.') - if trial_functions[0].space.codomain_complex==True: + if trial_functions[0].space.codomain_complex: raise TypeError('Trial space and Test space should be real. In the complex case, a SesquilinearForm should be called ') @@ -544,14 +546,14 @@ def __new__(cls, arguments, expr, **options): # Distinguish between trial and test functions trial_functions, test_functions = args - if trial_functions.space.codomain_complex != test_functions.space.codomain_complex: + if trial_functions[0].space.codomain_complex != test_functions[0].space.codomain_complex: raise TypeError('Trial space and Test space should both be complex.') - if trial_functions.space.codomain_complex==True: + if not trial_functions[0].space.codomain_complex: raise TypeError('Trial space and Test space should be complex. In the real case, a BilinearForm should be called ') # Check linearity with respect to trial functions - if not is_antilinear_expression(expr, trial_functions): + if not is_linear_expression(expr, trial_functions): msg = ' Expression is not linear w.r.t trial functions {}'\ .format(trial_functions) raise UnconsistentLinearExpressionError(msg) @@ -638,6 +640,7 @@ def __call__(self, trials, tests, **kwargs): #============================================================================== class Norm(Functional): is_norm = True + is_sesquilinear = False def __new__(cls, expr, domain, kind='l2', evaluate=True, **options): # # ... diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 5c8b6589..2bd41b9d 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -80,7 +80,7 @@ def get_logical_test_function(u): kind = space.kind dim = space.ldim logical_domain = space.domain.logical_domain - l_space = type(space)(space.name, logical_domain, kind=kind) + l_space = type(space)(space.name, logical_domain, kind=kind, codomain_complex=space.codomain_complex) el = l_space.element(u.name) return el @@ -895,7 +895,7 @@ def eval(cls, expr, domain, **options): """.""" from sympde.expr.evaluation import TerminalExpr, DomainExpression - from sympde.expr.expr import BilinearForm, LinearForm, BasicForm, Norm + from sympde.expr.expr import BilinearForm, LinearForm, BasicForm, Norm, SesquilinearForm from sympde.expr.expr import Integral types = (ScalarFunction, VectorFunction, DifferentialOperator, Trace, Integral) @@ -1221,6 +1221,12 @@ def eval(cls, expr, domain, **options): body = cls.eval(expr.expr, domain) return BilinearForm((trials, tests), body) + elif isinstance(expr, SesquilinearForm): + 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, domain) + return SesquilinearForm((trials, tests), body) + elif isinstance(expr, LinearForm): tests = [get_logical_test_function(a) for a in expr.test_functions] body = cls.eval(expr.expr, domain) From 9e4015aa2c09f44f0d8e8917a44c66d0873ef2e5 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Tue, 4 Jul 2023 09:00:33 +0200 Subject: [PATCH 5/9] make inner sustain 1D --- sympde/calculus/core.py | 9 ++++++++- sympde/expr/equation.py | 2 +- sympde/expr/expr.py | 7 +++++-- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index 2bb25312..e9a15754 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -513,9 +513,16 @@ def __new__(cls, arg1, arg2, **options): args_2 = [i for i in b if not i.is_commutative] c2 = [i for i in b if not i in args_2] + # args_2 = [i for args_2 in b if not i.is_commutative] + # c2 = [i for i in b if not i in args_2] + + c = Mul(*c1)*Mul(*c2) + if args_1==[] and args_2==[]: + return(c) + + a = reduce(mul, args_1) b = reduce(mul, args_2) - c = Mul(*c1)*Mul(*c2) if str(a) > str(b): a,b = b,a diff --git a/sympde/expr/equation.py b/sympde/expr/equation.py index 137043d5..66cb759c 100644 --- a/sympde/expr/equation.py +++ b/sympde/expr/equation.py @@ -11,7 +11,7 @@ from sympde.topology.space import VectorFunction, IndexedVectorFunction from sympde.topology import Boundary, NormalVector, TangentVector from sympde.topology import Trace, trace_0, trace_1 -from sympde.calculus import grad, dot +from sympde.calculus import grad, dot, inner from sympde.core.utils import random_string from .expr import BilinearForm, LinearForm, SesquilinearForm diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index b9a949c2..11952a25 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -350,7 +350,7 @@ def __new__(cls, arguments, expr, **options): args = _sanitize_arguments(arguments, is_linear=True) # Check linearity with respect to the given arguments - if not is_linear_expression(expr, args): + if not (is_linear_expression(expr, args) or is_antilinear_expression(expr, args)): msg = 'Expression is not linear w.r.t [{}]'.format(args) raise UnconsistentLinearExpressionError(msg) @@ -615,7 +615,7 @@ def is_hermitian(self): left, right = self.variables a1 = self(left, right) a2 = self(right, left) - self._is_hermitian = (a1 == a2) + self._is_hermitian = (conjugate(a1) == a2) return self._is_hermitian def __call__(self, trials, tests, **kwargs): @@ -849,6 +849,9 @@ def is_linear_expression(expr, args, integral=True, debug=True): if not( (newexpr-left_expr).expand() == 0 or newexpr.expand()==left_expr.expand()): # TODO use a warning or exception? if debug: + a=(newexpr-left_expr).expand() + b=newexpr.expand() + c=left_expr.expand() print('Failed to assert multiplication property') print('{} != {}'.format(newexpr, left_expr)) return False From d6621d890239cc486c453ae5b159dd2918db3b45 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Wed, 5 Jul 2023 09:02:23 +0200 Subject: [PATCH 6/9] Add is_complex property into ScalarFunction and VectorFunction + add sustainability od Abs function + Add is_sesquilinear flag to Functional --- sympde/expr/expr.py | 3 ++- sympde/topology/derivatives.py | 5 ++++- sympde/topology/space.py | 4 ++++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 11952a25..13f5d0a0 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -289,7 +289,8 @@ def _sympystr(self, printer): #============================================================================== class Functional(BasicForm): - is_functional = True + is_functional = True + is_sesquilinear = False def __new__(cls, expr, domain, evaluate=True, **options): diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index f6cdc688..caca41cf 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -15,7 +15,7 @@ from sympy import S from sympy import Indexed from sympy import diff -from sympy import log, conjugate +from sympy import log, conjugate, Abs from sympy import preorder_traversal from sympy import cacheit from sympy.core.compatibility import is_sequence @@ -165,6 +165,9 @@ def eval(cls, expr): elif isinstance(expr, conjugate): return conjugate(cls(expr.args[0], evaluate=True)) + elif isinstance(expr, Abs): + return Abs(cls(expr.args[0], evaluate=True)) + else: msg = '{expr} of type {type}'.format(expr=expr, type=type(expr)) raise NotImplementedError(msg) diff --git a/sympde/topology/space.py b/sympde/topology/space.py index 3f829b72..2b88570b 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -438,6 +438,7 @@ class ScalarFunction(Symbol): is_commutative = True _space = None _projection_of = None + is_complex = False def __new__(cls, space, name): if not isinstance(space, ScalarFunctionSpace): @@ -445,6 +446,7 @@ def __new__(cls, space, name): obj = Expr.__new__(cls) obj._space = space obj._name = name + obj.is_complex = space.codomain_complex return obj @property @@ -547,12 +549,14 @@ class VectorFunction(Symbol, IndexedBase): is_commutative = False _space = None _projection_of = None + is_complex = False def __new__(cls, space, name): if not isinstance(space, VectorFunctionSpace): raise ValueError('Expecting a VectorFunctionSpace') obj = Expr.__new__(cls) obj._space = space + obj.is_complex=space.codomain_complex obj._name = name return obj From a8a90e33e4c01d69e43fa95e1865b90732fa1409 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Wed, 5 Jul 2023 14:25:55 +0200 Subject: [PATCH 7/9] Replace definition of Norm + antilinearity + conjugate in Inner --- sympde/calculus/core.py | 17 ++++++++++++++--- sympde/expr/expr.py | 27 ++++++++++++++++----------- 2 files changed, 30 insertions(+), 14 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index e9a15754..c4b4e909 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -184,6 +184,17 @@ def is_zero(x): else: return x == 0 +def is_Function_complex(list): + res=False + for expr in list: + if isinstance(expr, (ScalarFunction, VectorFunction)): + res = res or expr.is_complex + else: + res = res or is_Function_complex(expr.args) + if res: + break + return res + #============================================================================== class BasicOperator(CalculusFunction): """ @@ -513,14 +524,14 @@ def __new__(cls, arg1, arg2, **options): args_2 = [i for i in b if not i.is_commutative] c2 = [i for i in b if not i in args_2] - # args_2 = [i for args_2 in b if not i.is_commutative] - # c2 = [i for i in b if not i in args_2] + if is_Function_complex(b): + args_2 = [conjugate(i) for i in args_2] + c2 = [conjugate(i) for i in c2] c = Mul(*c1)*Mul(*c2) if args_1==[] and args_2==[]: return(c) - a = reduce(mul, args_1) b = reduce(mul, args_2) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 13f5d0a0..ccabee2e 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -351,11 +351,16 @@ def __new__(cls, arguments, expr, **options): args = _sanitize_arguments(arguments, is_linear=True) # Check linearity with respect to the given arguments - if not (is_linear_expression(expr, args) or is_antilinear_expression(expr, args)): - msg = 'Expression is not linear w.r.t [{}]'.format(args) - raise UnconsistentLinearExpressionError(msg) + if args[0].is_complex: + if not is_antilinear_expression(expr, args): + msg = 'Complex case : Expression is not antilinear w.r.t [{}]'.format(args) + raise UnconsistentLinearExpressionError(msg) + else: + if not is_linear_expression(expr, args): + msg = 'Real case :Expression is not linear w.r.t [{}]'.format(args) + raise UnconsistentLinearExpressionError(msg) - # Create new object of type LinearForm + # Create new object of type LinearForm obj = Basic.__new__(cls, args, expr) # Compute 'domain' property (scalar or tuple) @@ -555,13 +560,13 @@ def __new__(cls, arguments, expr, **options): # Check linearity with respect to trial functions if not is_linear_expression(expr, trial_functions): - msg = ' Expression is not linear w.r.t trial functions {}'\ + msg = ' Real Case : Expression is not linear w.r.t trial functions {}'\ .format(trial_functions) raise UnconsistentLinearExpressionError(msg) # Check linearity with respect to test functions - if not is_linear_expression(expr, test_functions): - msg = ' Expression is not linear w.r.t test functions {}'\ + if not is_antilinear_expression(expr, test_functions): + msg = ' Complex Case : Expression is not linear w.r.t test functions {}'\ .format(test_functions) raise UnconsistentLinearExpressionError(msg) @@ -673,21 +678,21 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options): exponent = 2 if not is_vector: - expr = expr*expr + expr = Inner(expr,expr) else: if not( expr.shape[1] == 1 ): raise ValueError('Wrong expression for Matrix. must be a row') v = Tuple(*expr[:,0]) - expr = Dot(v, v) + expr = Inner(v, v) elif kind == 'h1'and evaluate : exponent = 2 if not is_vector: a = Grad(expr) - expr = Dot(a, a) + expr = Inner(a, a) else: if not( expr.shape[1] == 1 ): @@ -702,7 +707,7 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options): if not is_vector: a = Hessian(expr) - expr = Dot(a, a) + expr = Inner(a, a) else: raise NotImplementedError('TODO') From 602dca7510bf636a9cc3b6818107affe6dc8f358 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Mon, 10 Jul 2023 13:49:18 +0200 Subject: [PATCH 8/9] Put the conjugate intothe definition of Norm --- sympde/expr/expr.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index ccabee2e..e1c95d8b 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -678,36 +678,40 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options): exponent = 2 if not is_vector: - expr = Inner(expr,expr) + expr = expr * conjugate(expr) else: if not( expr.shape[1] == 1 ): raise ValueError('Wrong expression for Matrix. must be a row') - v = Tuple(*expr[:,0]) - expr = Inner(v, v) + args = expr[:, 0] + args_conj=[conjugate(i) for i in args] + v = Tuple(*args) + v_conj = Tuple(*args_conj) + expr = Dot(v, v_conj) elif kind == 'h1'and evaluate : exponent = 2 if not is_vector: a = Grad(expr) - expr = Inner(a, a) + expr = Dot(a, conjugate(a)) else: if not( expr.shape[1] == 1 ): raise ValueError('Wrong expression for Matrix. must be a row') - - v = Tuple(*expr[:,0]) - a = Grad(v) - expr = Inner(a, a) + args = expr[:, 0] + args_conj=[conjugate(i) for i in args] + v = Tuple(*args) + v_conj = Tuple(*args_conj) + expr = Dot(Grad(v), Grad(v_conj)) elif kind == 'h2'and evaluate : exponent = 2 if not is_vector: a = Hessian(expr) - expr = Inner(a, a) + expr = Dot(a, conjugate(a)) else: raise NotImplementedError('TODO') From eccc63aa046b6534609b699ab0503ea7463a1473 Mon Sep 17 00:00:00 2001 From: tomcaruso Date: Mon, 10 Jul 2023 14:03:22 +0200 Subject: [PATCH 9/9] Put conjugate of Inner in comment for now --- sympde/calculus/core.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index c4b4e909..c2629548 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -184,16 +184,16 @@ def is_zero(x): else: return x == 0 -def is_Function_complex(list): - res=False - for expr in list: - if isinstance(expr, (ScalarFunction, VectorFunction)): - res = res or expr.is_complex - else: - res = res or is_Function_complex(expr.args) - if res: - break - return res +# def is_Function_complex(list): +# res=False +# for expr in list: +# if isinstance(expr, (ScalarFunction, VectorFunction)): +# res = res or expr.is_complex +# else: +# res = res or is_Function_complex(expr.args) +# if res: +# break +# return res #============================================================================== class BasicOperator(CalculusFunction): @@ -524,9 +524,9 @@ def __new__(cls, arg1, arg2, **options): args_2 = [i for i in b if not i.is_commutative] c2 = [i for i in b if not i in args_2] - if is_Function_complex(b): - args_2 = [conjugate(i) for i in args_2] - c2 = [conjugate(i) for i in c2] + # if is_Function_complex(b): + # args_2 = [i if isinstance(i, conjugate) else conjugate(i) for i in args_2] + # c2 = [i if isinstance(i, conjugate) else conjugate(i) for i in c2] c = Mul(*c1)*Mul(*c2) if args_1==[] and args_2==[]: @@ -535,8 +535,8 @@ def __new__(cls, arg1, arg2, **options): a = reduce(mul, args_1) b = reduce(mul, args_2) - if str(a) > str(b): - a,b = b,a + # if str(a) > str(b): + # a,b = b,a obj = Basic.__new__(cls, a, b)