From a6e4faa746b1c17d7953b03f9dfd2f471a060fb1 Mon Sep 17 00:00:00 2001 From: ratnania Date: Mon, 19 Feb 2024 12:18:01 +0100 Subject: [PATCH 01/21] add Matrix object --- sympde/calculus/matrices.py | 58 +++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/sympde/calculus/matrices.py b/sympde/calculus/matrices.py index 16b827e1..08da3e2f 100644 --- a/sympde/calculus/matrices.py +++ b/sympde/calculus/matrices.py @@ -1,5 +1,6 @@ from sympy import Expr, S from sympy import Add, Mul, Pow +from sympy import Tuple from sympy import sympify from sympy.core.decorators import call_highest_priority from sympde.core.basic import _coeffs_registery, Basic @@ -289,6 +290,63 @@ def _sympystr(self, printer): sstr = printer.doprint return '{}[{}]'.format(sstr(self.args[0]),sstr(self.args[1])) +class Matrix(MatrixSymbolicExpr): + + def __new__(cls, *args, **options): + name = options.pop('name') + if name is None: + raise ValueError('Expecting a name keyword.') + + if not( len(args) == 1 ): + raise ValueError('Expecting one argument.') + + _args = args[0] + if not isinstance(_args, list): + raise TypeError('Expecting a list.') + + for _ in _args: + if not isinstance(_, list): + raise TypeError('args components must be a list.') + + sizes = [len(_) for _ in _args] + if not( len(set(sizes)) == 1): + raise ValueError('Wrong matrix size') + + # make _args a Tuple of Tuples + newargs = [] + for _ in _args: + a = Tuple(*_) + newargs.append(a) + _args = Tuple(*newargs) + + args = (_args, args[1:]) + + obj = Expr.__new__(cls, *args) + obj._name = name + + return obj + + @property + def name(self): + return self._name + + def __getitem__(self, key): + try: + i,j = key + return self.args[0][i][j] + except: + return MatrixElement(self, key) + + def _sympystr(self, printer): + sstr = printer.doprint + return '{}'.format(sstr(self.name)) + + def __hash__(self): + return hash((self.name, self.args)) + + + + Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = { "Mul": [lambda x: MatSymbolicMul(*x.args)], "Add": [lambda x: MatSymbolicAdd(*x.args)] From 9f8477353325ebac8c22d7476ef418392638f124 Mon Sep 17 00:00:00 2001 From: ratnania Date: Mon, 19 Feb 2024 13:08:53 +0100 Subject: [PATCH 02/21] add Vector class --- sympde/calculus/matrices.py | 43 +++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/sympde/calculus/matrices.py b/sympde/calculus/matrices.py index 08da3e2f..52a3d4c8 100644 --- a/sympde/calculus/matrices.py +++ b/sympde/calculus/matrices.py @@ -345,6 +345,49 @@ def __hash__(self): return hash((self.name, self.args)) +class Vector(MatrixSymbolicExpr): + + def __new__(cls, *args, **options): + name = options.pop('name') + if name is None: + raise ValueError('Expecting a name keyword.') + + if not( len(args) == 1 ): + raise ValueError('Expecting one argument.') + + _args = args[0] + if not isinstance(_args, list): + raise TypeError('Expecting a list.') + + # make _args a Tuple + _args = Tuple(*_args) + + args = (_args, args[1:]) + + obj = Expr.__new__(cls, *args) + obj._name = name + + return obj + + @property + def name(self): + return self._name + + def __getitem__(self, key): + try: + i = key + return self.args[0][i] + except: + # TODO ARA do we keep returning a MatrixElement? + return MatrixElement(self, key) + + def _sympystr(self, printer): + sstr = printer.doprint + return '{}'.format(sstr(self.name)) + + def __hash__(self): + return hash((self.name, self.args)) + Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = { From b9fe24e83528e8630909ba0bf988cc00c76b3f41 Mon Sep 17 00:00:00 2001 From: ratnania Date: Mon, 19 Feb 2024 14:17:15 +0100 Subject: [PATCH 03/21] mv matrices to core --- sympde/calculus/__init__.py | 1 - sympde/core/__init__.py | 1 + sympde/{calculus => core}/matrices.py | 0 sympde/expr/evaluation.py | 4 ++-- sympde/topology/mapping.py | 10 +++++----- 5 files changed, 8 insertions(+), 8 deletions(-) rename sympde/{calculus => core}/matrices.py (100%) diff --git a/sympde/calculus/__init__.py b/sympde/calculus/__init__.py index d468d7df..1a6a6448 100644 --- a/sympde/calculus/__init__.py +++ b/sympde/calculus/__init__.py @@ -1,3 +1,2 @@ from .core import * from .errors import * -from .matrices import * diff --git a/sympde/core/__init__.py b/sympde/core/__init__.py index 4bc8bacf..8d4e49aa 100644 --- a/sympde/core/__init__.py +++ b/sympde/core/__init__.py @@ -1,3 +1,4 @@ from .algebra import * from .basic import * from .utils import * +from .matrices import * diff --git a/sympde/calculus/matrices.py b/sympde/core/matrices.py similarity index 100% rename from sympde/calculus/matrices.py rename to sympde/core/matrices.py diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index e75189be..3a7f5c5c 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -20,12 +20,12 @@ Dot_2d, Inner_2d, Cross_2d, Dot_3d, Inner_3d, Cross_3d) from sympde.core.utils import random_string +from sympde.core.matrices import SymbolicDeterminant, Inverse, Transpose +from sympde.core.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace from sympde.calculus import jump, avg, minus, plus from sympde.calculus import Jump, is_zero from sympde.calculus.core import _generic_ops, _diff_ops -from sympde.calculus.matrices import SymbolicDeterminant, Inverse, Transpose -from sympde.calculus.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace from sympde.topology.basic import BasicDomain, Union, Interval from sympde.topology.basic import Boundary, Interface diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 5c8b6589..c5fa7908 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -19,12 +19,12 @@ from sympde.core.basic import BasicMapping from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery +from sympde.core.matrices import MatrixSymbolicExpr, MatrixElement, SymbolicTrace, Inverse +from sympde.core.matrices import SymbolicDeterminant, Transpose from sympde.calculus.core import PlusInterfaceOperator, MinusInterfaceOperator from sympde.calculus.core import grad, div, curl, laplace #, hessian from sympde.calculus.core import dot, inner, outer, _diff_ops from sympde.calculus.core import has, DiffOperator -from sympde.calculus.matrices import MatrixSymbolicExpr, MatrixElement, SymbolicTrace, Inverse -from sympde.calculus.matrices import SymbolicDeterminant, Transpose from .basic import BasicDomain, Union, InteriorDomain from .basic import Boundary, Connectivity, Interface @@ -976,7 +976,7 @@ def eval(cls, expr, domain, **options): elif isinstance(expr, Transpose): arg = cls(expr.arg, domain) return Transpose(arg) - + elif isinstance(expr, grad): arg = expr.args[0] if isinstance(mapping, InterfaceMapping): @@ -1117,7 +1117,7 @@ def eval(cls, expr, domain, **options): elif dim == 3: lgrad_arg = LogicalGrad_3d(arg) - + grad_arg = Covariant(mapping, lgrad_arg) expr = grad_arg[0] return expr @@ -1243,7 +1243,7 @@ def eval(cls, expr, domain, **options): domain = domain.logical_domain det = TerminalExpr(sqrt((J.T*J).det()), domain=domain) return DomainExpression(domain, ImmutableDenseMatrix([[newexpr*det]])) - + elif isinstance(expr, Function): args = [cls.eval(a, domain) for a in expr.args] return type(expr)(*args) From 572011d313afbf32f8c422a51514ce2dd0b10b0d Mon Sep 17 00:00:00 2001 From: ratnania Date: Mon, 19 Feb 2024 19:40:26 +0100 Subject: [PATCH 04/21] TerminalExpr for sympde objects: Matrix and Vector --- sympde/expr/evaluation.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 3a7f5c5c..3bc3dbd7 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -22,6 +22,8 @@ from sympde.core.utils import random_string from sympde.core.matrices import SymbolicDeterminant, Inverse, Transpose from sympde.core.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace +from sympde.core.matrices import Matrix as SympdeMatrix +from sympde.core.matrices import Vector as SympdeVector from sympde.calculus import jump, avg, minus, plus from sympde.calculus import Jump, is_zero @@ -614,6 +616,12 @@ def eval(cls, expr, domain): elif isinstance(expr, Inverse): return cls.eval(expr.arg, domain=domain).inv() + elif isinstance(expr, SympdeMatrix): + return ImmutableDenseMatrix(expr.args[0]) + + elif isinstance(expr, SympdeVector): + return ImmutableDenseMatrix([expr.args[0]]) + elif isinstance(expr, ScalarFunction): return expr From 3a658bcd4e12578da6137ddc3f67a9cca4e6ec27 Mon Sep 17 00:00:00 2001 From: ratnania Date: Wed, 21 Feb 2024 10:59:49 +0100 Subject: [PATCH 05/21] add __eq__ to vector and matrix --- sympde/core/matrices.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index 52a3d4c8..b530b17f 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -344,6 +344,12 @@ def _sympystr(self, printer): def __hash__(self): return hash((self.name, self.args)) + def __eq__(self, a): + if isinstance(a, Matrix): + eq = self.name == a.name + return eq + return False + class Vector(MatrixSymbolicExpr): @@ -388,6 +394,12 @@ def _sympystr(self, printer): def __hash__(self): return hash((self.name, self.args)) + def __eq__(self, a): + if isinstance(a, Matrix): + eq = self.name == a.name + return eq + return False + Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = { From d02fd51999841bd7539a2fdbb558398327f91dcf Mon Sep 17 00:00:00 2001 From: ratnania Date: Wed, 21 Feb 2024 11:00:10 +0100 Subject: [PATCH 06/21] add new tests for matrices and vectors --- sympde/expr/tests/test_expr.py | 275 ++++++++++++++++++++++++++++----- 1 file changed, 237 insertions(+), 38 deletions(-) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index efca3a98..62aca72d 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -10,9 +10,13 @@ from sympy import ImmutableDenseMatrix as Matrix from sympde.core import Constant -from sympde.calculus import grad, dot, inner, rot, div +from sympde.calculus import dot, inner, outer, cross +from sympde.calculus import grad, rot, div, curl from sympde.calculus import laplace, bracket, convect from sympde.calculus import jump, avg, Dn, minus, plus +from sympde.core import Matrix as SympdeMatrix +from sympde.core import Vector +from sympde.core import Transpose from sympde.topology import dx1, dx2, dx3 from sympde.topology import dx, dy, dz @@ -33,6 +37,8 @@ from sympde.expr.expr import Functional, Norm from sympde.expr.expr import linearize from sympde.expr.evaluation import TerminalExpr +from sympde.expr.expr import is_linear_expression + #============================================================================== def test_linear_expr_2d_1(): @@ -213,7 +219,7 @@ def test_linear_form_2d_2(): int_0 = lambda expr: integral(domain , expr) int_1 = lambda expr: integral(B1, expr) - g = Matrix((x,y)) + g = Vector([x,y], name='g') l = LinearForm(v, int_0(dot(g, v))) assert(l.domain == domain.interior) @@ -221,7 +227,7 @@ def test_linear_form_2d_2(): # ... # ... - g = Matrix((x,y)) + g = Vector([x,y], name='g') l1 = LinearForm(v1, int_0(dot(g, v1))) l = LinearForm(v, l1(v)) @@ -230,8 +236,8 @@ def test_linear_form_2d_2(): # ... # ... - g1 = Matrix((x,0)) - g2 = Matrix((0,y)) + g1 = Vector([x,0], name='g1') + g2 = Vector([0,y], name='g2') l1 = LinearForm(v1, int_0(dot(v1, g1))) l2 = LinearForm(v2, int_0(dot(v2, g2))) @@ -396,14 +402,14 @@ def test_terminal_expr_linear_2d_1(): # ... # ... - g = Matrix((x**2, y**2)) + g = Vector([x**2, y**2], name='g') l = LinearForm(v, int_1(v*dot(g, nn))) print(TerminalExpr(l, domain)) print('') # ... # ... - g = Matrix((x**2, y**2)) + g = Vector([x**2, y**2], name='g') l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v)) print(TerminalExpr(l, domain)) print('') @@ -417,7 +423,7 @@ def test_terminal_expr_linear_2d_1(): # ... # ... - g = Matrix((x,y)) + g = Vector([x,y], name='g') l1 = LinearForm(v1, int_0(x*y*v1)) l2 = LinearForm(v2, int_0(dot(grad(v2), g))) @@ -435,7 +441,7 @@ def test_terminal_expr_linear_2d_1(): # ... # ... - g = Matrix((x**2, y**2)) + g = Vector([x**2, y**2], name='g') l1 = LinearForm(v1, int_0(x*y*v1)) l2 = LinearForm(v1, int_0(v1)) l3 = LinearForm(v, int_1(v*dot(g, nn))) @@ -464,14 +470,14 @@ def test_terminal_expr_linear_2d_2(): int_0 = lambda expr: integral(domain , expr) int_1 = lambda expr: integral(B1, expr) - g = Matrix((x,y)) + g = Vector([x, y], name='g') l = LinearForm(v, int_0(dot(g, v))) print(TerminalExpr(l, domain)) print('') # ... # ... - g = Matrix((x,y)) + g = Vector([x, y], name='g') l = LinearForm(v, int_0(dot(g, v) + div(v))) print(TerminalExpr(l, domain)) print('') @@ -585,7 +591,7 @@ def test_terminal_expr_linear_2d_5(): V = ScalarFunctionSpace('V', domain) - B_neumann = [domain.get_boundary(axis=0, ext=-1), + B_neumann = [domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=1, ext=-1)] if len(B_neumann) == 1: @@ -1178,9 +1184,7 @@ def test_stabilization_2d_1(): kappa = Constant('kappa', is_real=True) mu = Constant('mu' , is_real=True) - b1 = 1. - b2 = 0. - b = Matrix((b1, b2)) + b = Vector([1., 0.], name='b') # right hand side f = x*y @@ -1393,25 +1397,26 @@ def test_norm_2d_1(): print('') # ... -#============================================================================== -def test_norm_2d_2(): - - domain = Domain('Omega', dim=2) - x,y = domain.coordinates - - V = VectorFunctionSpace('V', domain) - F = element_of(V, 'F') - - # ... - f = Matrix((sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y))) - expr = Matrix([F[0]-f[0], F[1]-f[1]]) - l2_norm_u = Norm(expr, domain, kind='l2') - h1_norm_u = Norm(expr, domain, kind='h1') - - print('> l2 norm = ', TerminalExpr(l2_norm_u, domain)) - print('> h1 norm = ', TerminalExpr(h1_norm_u, domain)) - print('') - # ... +##============================================================================== +## TODO ARA to be fixed +#def test_norm_2d_2(): +# +# domain = Domain('Omega', dim=2) +# x,y = domain.coordinates +# +# V = VectorFunctionSpace('V', domain) +# F = element_of(V, 'F') +# +# # ... +# f = Vector([sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y)], name='f') +# expr = Vector([F[0]-f[0], F[1]-f[1]], name='expr') +# l2_norm_u = Norm(expr, domain, kind='l2') +# h1_norm_u = Norm(expr, domain, kind='h1') +# +# print('> l2 norm = ', TerminalExpr(l2_norm_u, domain)) +# print('> h1 norm = ', TerminalExpr(h1_norm_u, domain)) +# print('') +# # ... #============================================================================== # this test checks some properties of bilinear forms @@ -1502,16 +1507,16 @@ def test_linearity_linear_form_2d_1(): _ = LinearForm(v, int_0(x * y * v)) _ = LinearForm(v, int_0(x * y * v + v)) - g = Matrix((x**2, y**2)) + g = Vector([x**2, y**2], name='g') _ = LinearForm(v, int_0(v * dot(g, nn))) - g = Matrix((x**2, y**2)) + g = Vector([x**2, y**2], name='g') _ = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v)) l1 = LinearForm(v1, int_0(x * y * v1)) _ = LinearForm(v, l1(v)) - g = Matrix((x,y)) + g = Vector([x, y], name='g') l1 = LinearForm(v1, int_0(x*y*v1)) l2 = LinearForm(v2, int_0(dot(grad(v2), g))) _ = LinearForm(v, l1(v) + l2(v)) @@ -1520,7 +1525,7 @@ def test_linearity_linear_form_2d_1(): l2 = LinearForm(v1, int_0(v1)) _ = LinearForm(v, l1(v) + kappa*l2(v)) - g = Matrix((x**2, y**2)) + g = Vector([x**2, y**2], name='g') l1 = LinearForm(v1, int_0(x*y*v1)) l2 = LinearForm(v1, int_0(v1)) l3 = LinearForm(v, int_1(v*dot(g, nn))) @@ -1891,6 +1896,200 @@ def test_interface_integral_4(): assert expr[0].expr[1,1] == u2[1]*v2[1] # ... + + +#============================================================================== +def test_matrices_vectors(): + + domain = Domain('Omega', dim=3) + x,y,z = domain.coordinates + + kappa = Constant('kappa', is_real=True) + mu = Constant('mu', is_real=True) + theta = Constant('theta', is_real=True) + a1 = Constant('a1', is_real=True) + a2 = Constant('a2', is_real=True) + b1 = Constant('b1', is_real=True) + b2 = Constant('b2', is_real=True) + + V = ScalarFunctionSpace('V', domain) + W = VectorFunctionSpace('W', domain) + + u,v = [element_of(V, name=i) for i in ['u', 'v']] + F,G = [element_of(W, name=i) for i in ['F', 'G']] + + #A = SympdeMatrix([[cos(theta), 0], [0, sin(theta)]], name='A') + A = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A') + b = Vector([1, 2, 3], name='b') + + ######################################################################## + # SCALAR-FUNCTION CASE + ######################################################################## + # ... + f = lambda a: A*grad(u) + + assert is_linear_expression(f(u), (u,)) + # ... + + # ... + f = lambda u: b*u + + assert is_linear_expression(f(u), (u,)) + # ... + + # ... + f = lambda u: dot(b,grad(u)) + + assert is_linear_expression(f(u), (u,)) + # ... + + # ... + f = lambda u,v: dot(A*grad(u), grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(grad(u), A*grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(A*grad(u), A*grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(A.T*grad(u), grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(grad(u), A.T*grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(A.T*grad(u), A.T*grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(A*grad(u), grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(b,grad(u))*v + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(A*grad(u), grad(v)) + dot(b,grad(u))*v + dot(b,grad(v))*u + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: cross(b*u, b*v) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: cross(b*u, grad(v)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: cross(grad(u), b*v) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda u,v: dot(cross(grad(u), b), cross(grad(v), b)) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... + + # ... + f = lambda F,G: dot(curl(cross(b, F)), cross(b, G)) + + assert is_linear_expression(f(F,G), (F,)) + assert is_linear_expression(f(F,G), (G,)) + # ... + + ## ... + #f = lambda u,v: 0 + # + #assert is_linear_expression(f(u,v), (u,)) + #assert is_linear_expression(f(u,v), (v,)) + ## ... + # + ## ... + #f = lambda u,v: 0 + # + #assert is_linear_expression(f(u,v), (u,)) + #assert is_linear_expression(f(u,v), (v,)) + ## ... + + ######################################################################## + # VECTOR-FUNCTION CASE + ######################################################################## + # ... + f = lambda F,G: cross(A*F, G) + + assert is_linear_expression(f(F,G), (F,)) + assert is_linear_expression(f(F,G), (G,)) + # ... + + # ... + f = lambda F,G: cross(F, A*G) + + assert is_linear_expression(f(F,G), (F,)) + assert is_linear_expression(f(F,G), (G,)) + # ... + + # ... + f = lambda F,G: cross(A*F, A*G) + + assert is_linear_expression(f(F,G), (F,)) + assert is_linear_expression(f(F,G), (G,)) + # ... + + # ... + I = SympdeMatrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], name='I') + epsilon = lambda w: 0.5*(grad(w) + Transpose(grad(w))) + sigma = lambda w: kappa * div(w) * I + 2 * mu * epsilon(w) + f = lambda u,v: inner(sigma(u), epsilon(v)) + + assert is_linear_expression(f(F,G), (F,)) + assert is_linear_expression(f(F,G), (G,)) + # ... + + #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== From 760fdd17a48bf8ade4eaf9e9918534f5bee45384 Mon Sep 17 00:00:00 2001 From: ratnania Date: Wed, 21 Feb 2024 11:38:30 +0100 Subject: [PATCH 07/21] add to_tympy method that converts Matrix/Vector to ImmutableDenseMatrix from sympy --- sympde/core/matrices.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index b530b17f..2f6ca676 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -5,6 +5,9 @@ from sympy.core.decorators import call_highest_priority from sympde.core.basic import _coeffs_registery, Basic +from sympy import ImmutableDenseMatrix + + class MatrixSymbolicExpr(Expr): is_commutative = False _op_priority = 20.0 @@ -350,6 +353,9 @@ def __eq__(self, a): return eq return False + def to_sympy(self): + return ImmutableDenseMatrix(self.args[0]) + class Vector(MatrixSymbolicExpr): @@ -400,6 +406,9 @@ def __eq__(self, a): return eq return False + def to_sympy(self): + _args = [[_] for _ in self.args[0]] + return ImmutableDenseMatrix(_args) Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = { From e51563c61c1f78c194a6bf30c9dd141dc31a8f94 Mon Sep 17 00:00:00 2001 From: ratnania Date: Wed, 21 Feb 2024 11:39:21 +0100 Subject: [PATCH 08/21] fixing test_norm_2d_2 test, for Norm --- sympde/expr/tests/test_expr.py | 39 +++++++++++++++++----------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 62aca72d..2c09f00a 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -1397,26 +1397,25 @@ def test_norm_2d_1(): print('') # ... -##============================================================================== -## TODO ARA to be fixed -#def test_norm_2d_2(): -# -# domain = Domain('Omega', dim=2) -# x,y = domain.coordinates -# -# V = VectorFunctionSpace('V', domain) -# F = element_of(V, 'F') -# -# # ... -# f = Vector([sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y)], name='f') -# expr = Vector([F[0]-f[0], F[1]-f[1]], name='expr') -# l2_norm_u = Norm(expr, domain, kind='l2') -# h1_norm_u = Norm(expr, domain, kind='h1') -# -# print('> l2 norm = ', TerminalExpr(l2_norm_u, domain)) -# print('> h1 norm = ', TerminalExpr(h1_norm_u, domain)) -# print('') -# # ... +#============================================================================== +def test_norm_2d_2(): + + domain = Domain('Omega', dim=2) + x,y = domain.coordinates + + V = VectorFunctionSpace('V', domain) + F = element_of(V, 'F') + + # ... + f = Vector([sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y)], name='f') + expr = Vector([F[0]-f[0], F[1]-f[1]], name='expr') + l2_norm_u = Norm(expr, domain, kind='l2') + h1_norm_u = Norm(expr, domain, kind='h1') + + print('> l2 norm = ', TerminalExpr(l2_norm_u, domain)) + print('> h1 norm = ', TerminalExpr(h1_norm_u, domain)) + print('') + # ... #============================================================================== # this test checks some properties of bilinear forms From d9bcba7f5736e78dda164ac84e5756c3f8e542c3 Mon Sep 17 00:00:00 2001 From: ratnania Date: Wed, 21 Feb 2024 11:41:05 +0100 Subject: [PATCH 09/21] convert Norm expression Matrix/Vector to sympy --- sympde/expr/expr.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index bd87dfba..c866fc8a 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -22,6 +22,8 @@ from sympde.topology.space import ScalarFunction from sympde.topology.space import VectorFunction from sympde.topology.space import Trace, trace_0, trace_1 +from sympde.core.matrices import Matrix as SympdeMatrix +from sympde.core.matrices import Vector from .errors import UnconsistentLinearExpressionError from .basic import BasicForm @@ -538,9 +540,12 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options): # ... # ... - is_vector = isinstance(expr, (Matrix, ImmutableDenseMatrix, Tuple, list, tuple)) + is_vector = isinstance(expr, (Vector, SympdeMatrix, Matrix, ImmutableDenseMatrix, Tuple, list, tuple)) if is_vector: - expr = ImmutableDenseMatrix(expr) + if isinstance(expr, (Vector, SympdeMatrix)): + expr = expr.to_sympy() + else: + expr = ImmutableDenseMatrix(expr) # ... # ... @@ -610,9 +615,12 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options): # ... # ... - is_vector = isinstance(expr, (Matrix, ImmutableDenseMatrix, Tuple, list, tuple)) + is_vector = isinstance(expr, (Vector, SympdeMatrix, Matrix, ImmutableDenseMatrix, Tuple, list, tuple)) if is_vector: - expr = ImmutableDenseMatrix(expr) + if isinstance(expr, (Vector, SympdeMatrix)): + expr = expr.to_sympy() + else: + expr = ImmutableDenseMatrix(expr) # ... # ... From d1928b40d56ef7a8dac6b2887dcf4fdf7c7658ba Mon Sep 17 00:00:00 2001 From: ratnania Date: Wed, 21 Feb 2024 11:46:07 +0100 Subject: [PATCH 10/21] fix evaluation of Vector in TerminalExpr --- sympde/expr/evaluation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 3bc3dbd7..0d2ae4fa 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -617,10 +617,11 @@ def eval(cls, expr, domain): return cls.eval(expr.arg, domain=domain).inv() elif isinstance(expr, SympdeMatrix): - return ImmutableDenseMatrix(expr.args[0]) + return expr.to_sympy() elif isinstance(expr, SympdeVector): - return ImmutableDenseMatrix([expr.args[0]]) +# return ImmutableDenseMatrix([expr.args[0]]) + return expr.to_sympy() elif isinstance(expr, ScalarFunction): return expr From bdfe6eb0ad8d3cb96f976d4049907b51c5ec4509 Mon Sep 17 00:00:00 2001 From: ratnania Date: Thu, 22 Feb 2024 17:53:26 +0100 Subject: [PATCH 11/21] update mapping and derivatives to use Vector as input instead of Tuple. Update all examples and remove the use of Tuple as an input --- sympde/expr/tests/test_equation_2d.py | 13 ++++++------- sympde/expr/tests/test_expr.py | 22 +++++++++++----------- sympde/expr/tests/test_tensor_2d.py | 4 ++-- sympde/exterior/tests/test_datatype.py | 1 - sympde/exterior/tests/test_exterior.py | 1 - sympde/exterior/tests/test_inference.py | 1 - sympde/topology/derivatives.py | 6 ++++++ sympde/topology/mapping.py | 12 ++++++++++-- sympde/topology/tests/test_derivatives.py | 21 +++++++++++---------- sympde/topology/tests/test_gallery.py | 1 - sympde/topology/tests/test_mapping.py | 7 ++++--- 11 files changed, 50 insertions(+), 39 deletions(-) diff --git a/sympde/expr/tests/test_equation_2d.py b/sympde/expr/tests/test_equation_2d.py index cb446f2f..5d66e10e 100644 --- a/sympde/expr/tests/test_equation_2d.py +++ b/sympde/expr/tests/test_equation_2d.py @@ -2,7 +2,6 @@ #import pytest -from sympy.core.containers import Tuple from sympy import pi, cos, sin from sympy import ImmutableDenseMatrix as Matrix @@ -48,7 +47,7 @@ def test_equation_2d_1(): int_0 = lambda expr: integral(domain , expr) int_1 = lambda expr: integral(B1, expr) int_2 = lambda expr: integral(B2, expr) - + # ... bilinear/linear forms expr = dot(grad(v), grad(u)) a1 = BilinearForm((v,u), int_0(expr)) @@ -169,7 +168,7 @@ def test_equation_2d_2(): alpha = Constant('alpha', real=True) int_0 = lambda expr: integral(domain , expr) - + s = BilinearForm((tau,sigma), int_0(dot(grad(tau), grad(sigma)))) m = BilinearForm((tau,sigma), int_0(tau*sigma)) b1 = BilinearForm((tau,dw), int_0(bracket(pn, dw) * tau)) @@ -202,7 +201,7 @@ def test_equation_2d_3(): x,y = domain.coordinates B1 = Boundary(r'\Gamma_1', domain) - + int_0 = lambda expr: integral(domain , expr) int_1 = lambda expr: integral(B1, expr) @@ -235,7 +234,7 @@ def test_equation_2d_4(): x,y = domain.coordinates B1 = Boundary(r'\Gamma_1', domain) - + int_0 = lambda expr: integral(domain , expr) int_1 = lambda expr: integral(B1, expr) @@ -283,7 +282,7 @@ def test_equation_2d_5(): p,q = [element_of(V, name=i) for i in ['p', 'q']] int_0 = lambda expr: integral(domain , expr) - + a0 = BilinearForm((v,u), int_0(inner(grad(v), grad(u)))) print(' a0 done.') a1 = BilinearForm((q,p), int_0(p*q)) @@ -335,7 +334,7 @@ def test_equation_2d_6(): u,v = [element_of(V, name=i) for i in ['u', 'v']] int_0 = lambda expr: integral(domain , expr) - + # ... expr = kappa * dot(grad(u), grad(v)) + dot(b, grad(u)) * v a = BilinearForm((v,u), int_0(expr)) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 2c09f00a..7dca172f 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -4,7 +4,6 @@ import pytest -from sympy.core.containers import Tuple from sympy import Function from sympy import pi, cos, sin, exp from sympy import ImmutableDenseMatrix as Matrix @@ -17,6 +16,7 @@ from sympde.core import Matrix as SympdeMatrix from sympde.core import Vector from sympde.core import Transpose +from sympde.core import Vector from sympde.topology import dx1, dx2, dx3 from sympde.topology import dx, dy, dz @@ -88,7 +88,7 @@ def test_linear_expr_2d_2(): u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']] v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']] - g = Tuple(x,y) + g = Vector([x,y], name='g') l = LinearExpr(v, dot(g, v)) print(l) print(l.expr) @@ -99,8 +99,8 @@ def test_linear_expr_2d_2(): # ... # ... - g1 = Tuple(x,0) - g2 = Tuple(0,y) + g1 = Vector([x,0], name='g1') + g2 = Vector([0,y], name='g2') l = LinearExpr((v1,v2), dot(g1, v1) + dot(g2, v2)) print(l) print(l.expr) @@ -138,7 +138,7 @@ def test_linear_form_2d_1(): # ... # ... - g = Tuple(x**2, y**2) + g = Vector([x**2, y**2], name='g') l = LinearForm(v, int_1(v*dot(g, nn))) print(l) @@ -148,7 +148,7 @@ def test_linear_form_2d_1(): # ... # ... - g = Tuple(x**2, y**2) + g = Vector([x**2, y**2], name='g') l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v)) assert(len(l.domain.args) == 2) @@ -168,7 +168,7 @@ def test_linear_form_2d_1(): # ... # ... - g = Tuple(x,y) + g = Vector([x,y], name='g') l1 = LinearForm(u1, int_0(x*y*u1)) l2 = LinearForm(u2, int_0(dot(grad(u2), g))) @@ -485,14 +485,14 @@ def test_terminal_expr_linear_2d_2(): # TODO # # ... -# g = Tuple(x**2, y**2) +# g = Vector([x**2,y**2], name='g') # l = LinearForm(v, v*trace_1(g, B1)) # print(TerminalExpr(l)) # print('') # # ... # # # ... -# g = Tuple(x**2, y**2) +# g = Vector([x**2,y**2], name='g') # l = LinearForm(v, v*trace_1(g, B1) + x*y*v) # print(TerminalExpr(l)) # print('') @@ -506,7 +506,7 @@ def test_terminal_expr_linear_2d_2(): # # ... # # # ... -# g = Tuple(x,y) +# g = Vector([x,y], name='g') # l1 = LinearForm(v1, x*y*v1) # l2 = LinearForm(v2, dot(grad(v2), g)) # @@ -524,7 +524,7 @@ def test_terminal_expr_linear_2d_2(): # # ... # # # ... -# g = Tuple(x**2, y**2) +# g = Vector([x**2,y**2], name='g') # l1 = LinearForm(v1, x*y*v1) # l2 = LinearForm(v1, v1) # l3 = LinearForm(v, v*trace_1(g, B1)) diff --git a/sympde/expr/tests/test_tensor_2d.py b/sympde/expr/tests/test_tensor_2d.py index 0fc9e31d..33764286 100644 --- a/sympde/expr/tests/test_tensor_2d.py +++ b/sympde/expr/tests/test_tensor_2d.py @@ -2,9 +2,9 @@ # TODO: - add assert to every test -from sympy.core.containers import Tuple from sympde.core import Constant +from sympde.core import Vector from sympde.calculus import grad, dot, curl, div #from sympde.calculus import laplace #from sympde.topology import dx @@ -112,7 +112,7 @@ def test_tensorize_2d_3(): bx = Constant('bx') by = Constant('by') - b = Tuple(bx, by) + b = Vector([bx,by], name='b') expr = integral(domain, dot(b, grad(v)) * dot(b, grad(u))) a = BilinearForm((u,v), expr) diff --git a/sympde/exterior/tests/test_datatype.py b/sympde/exterior/tests/test_datatype.py index 9bea28d1..6a904464 100644 --- a/sympde/exterior/tests/test_datatype.py +++ b/sympde/exterior/tests/test_datatype.py @@ -2,7 +2,6 @@ from sympy import symbols -from sympy import Tuple from sympy import Matrix from sympy import srepr from sympy import Symbol diff --git a/sympde/exterior/tests/test_exterior.py b/sympde/exterior/tests/test_exterior.py index 07d28c88..caa882df 100644 --- a/sympde/exterior/tests/test_exterior.py +++ b/sympde/exterior/tests/test_exterior.py @@ -2,7 +2,6 @@ from sympy import symbols -from sympy import Tuple from sympy import Matrix from sympy import srepr from sympy import Symbol diff --git a/sympde/exterior/tests/test_inference.py b/sympde/exterior/tests/test_inference.py index 850cdaaa..ca12b2e3 100644 --- a/sympde/exterior/tests/test_inference.py +++ b/sympde/exterior/tests/test_inference.py @@ -2,7 +2,6 @@ from sympy import symbols -from sympy import Tuple from sympy import Matrix from sympy import srepr from sympy import Symbol diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index f5d3537a..2bc3f8a0 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -23,6 +23,8 @@ from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery from sympde.core.basic import BasicMapping +from sympde.core.matrices import Vector +from sympde.core.matrices import Matrix as SympdeMatrix from sympde.core.algebra import LinearOperator from sympde.calculus.core import minus, plus from sympde.calculus.core import has @@ -82,6 +84,10 @@ def eval(cls, expr): args = [cls(i, evaluate=True) for i in expr] args = Tuple(*args) return Matrix([args]) + elif isinstance(expr, Vector): + args = [[cls(_, evaluate=True)] for _ in expr.args[0]] + return ImmutableDenseMatrix(args) + elif isinstance(expr, (Matrix, ImmutableDenseMatrix)): newexpr = Matrix.zeros(*expr.shape) for i in range(expr.shape[0]): diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index c5fa7908..75ad3867 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -21,6 +21,8 @@ from sympde.core.basic import _coeffs_registery from sympde.core.matrices import MatrixSymbolicExpr, MatrixElement, SymbolicTrace, Inverse from sympde.core.matrices import SymbolicDeterminant, Transpose +from sympde.core.matrices import Vector +from sympde.core.matrices import Matrix as SympdeMatrix from sympde.calculus.core import PlusInterfaceOperator, MinusInterfaceOperator from sympde.calculus.core import grad, div, curl, laplace #, hessian from sympde.calculus.core import dot, inner, outer, _diff_ops @@ -797,7 +799,10 @@ def eval(cls, F, v): the covariant transformation """ - if not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)): + if isinstance(v, (Vector, SympdeMatrix)): + v = v.to_sympy() + + elif not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)): raise TypeError('> Expecting a tuple, list, Tuple, Matrix') assert F.pdim == F.ldim @@ -849,7 +854,10 @@ def eval(cls, F, v): if not isinstance(F, Mapping): raise TypeError('> Expecting a Mapping') - if not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)): + if isinstance(v, (Vector, SympdeMatrix)): + v = v.to_sympy() + + elif not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)): raise TypeError('> Expecting a tuple, list, Tuple, Matrix') M = Jacobian(F) diff --git a/sympde/topology/tests/test_derivatives.py b/sympde/topology/tests/test_derivatives.py index 3694adfe..b7cf31d0 100644 --- a/sympde/topology/tests/test_derivatives.py +++ b/sympde/topology/tests/test_derivatives.py @@ -1,11 +1,11 @@ # coding: utf-8 from sympy import symbols -from sympy import Tuple from sympy import Matrix from sympy import srepr from sympde.core import Constant +from sympde.core import Vector from sympde.calculus import grad, dot, inner from sympde.topology import Domain, element_of from sympde.topology import get_index_derivatives_atom @@ -25,6 +25,7 @@ def indices_as_str(a): # ... +# TODO ARA add more tests for sympde/Vector and Matrix def test_partial_derivatives_1(): print('============ test_partial_derivatives_1 ==============') @@ -39,7 +40,7 @@ def test_partial_derivatives_1(): V = ScalarFunctionSpace('V', mapped_domain) F,u,v,w = [element_of(V, name=i) for i in ['F', 'u', 'v', 'w']] - uvw = Tuple(u,v,w) + uvw = Vector([u,v,w], name='uvw') alpha = Constant('alpha') beta = Constant('beta') @@ -54,14 +55,14 @@ def test_partial_derivatives_1(): assert(dz(y**2) == 0) assert(dx(x*F) == F + x*dx(F)) - assert(dx(uvw) == Matrix([[dx(u), dx(v), dx(w)]])) - assert(dx(uvw) + dy(uvw) == Matrix([[dx(u) + dy(u), - dx(v) + dy(v), - dx(w) + dy(w)]])) - - expected = Matrix([[alpha*dx(u) + beta*dy(u), - alpha*dx(v) + beta*dy(v), - alpha*dx(w) + beta*dy(w)]]) + assert(dx(uvw) == Matrix([[dx(u)], [dx(v)], [dx(w)]])) + assert(dx(uvw) + dy(uvw) == Matrix([[dx(u) + dy(u)], + [dx(v) + dy(v)], + [dx(w) + dy(w)]])) + + expected = Matrix([[alpha*dx(u) + beta*dy(u)], + [alpha*dx(v) + beta*dy(v)], + [alpha*dx(w) + beta*dy(w)]]) assert(alpha * dx(uvw) + beta * dy(uvw) == expected) # ... diff --git a/sympde/topology/tests/test_gallery.py b/sympde/topology/tests/test_gallery.py index 6222c9ce..bf2875c1 100644 --- a/sympde/topology/tests/test_gallery.py +++ b/sympde/topology/tests/test_gallery.py @@ -1,7 +1,6 @@ # coding: utf-8 from sympy.abc import x,y,z -from sympy import Tuple from sympy import symbols x1, x2, x3 = symbols('x1, x2, x3') diff --git a/sympde/topology/tests/test_mapping.py b/sympde/topology/tests/test_mapping.py index 84ede14e..30e09df5 100644 --- a/sympde/topology/tests/test_mapping.py +++ b/sympde/topology/tests/test_mapping.py @@ -5,12 +5,13 @@ from sympy.tensor import IndexedBase from sympy import symbols, simplify +from sympde.core.matrices import Vector from sympde.topology import Mapping, MappedDomain from sympde.topology import dx, dy, dz from sympde.topology import dx1, dx2, dx3 from sympde.topology import Domain - from sympde.topology.mapping import Jacobian, Covariant, Contravariant + # ... def test_mapping_1d(): print('============ test_mapping_1d ==============') @@ -41,7 +42,7 @@ def test_mapping_2d(): F = Mapping('F', dim=dim) a,b = symbols('a b') - ab = Tuple(a, b) + ab = Vector([a, b], name='ab') assert(F.name == 'F') @@ -81,7 +82,7 @@ def test_mapping_3d(): F = Mapping('F', dim=dim) a,b,c = symbols('a b c') - abc = Tuple(a, b, c) + abc = Vector([a, b, c], name='abc') assert(F.name == 'F') From 369e295b60a7fa5ca8e39f64091dd39fb56d849e Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 10:38:08 +0100 Subject: [PATCH 12/21] we should only name a variable as _ if it is a local variable which is not used later --- sympde/topology/derivatives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index 2bc3f8a0..14169e7d 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -85,7 +85,7 @@ def eval(cls, expr): args = Tuple(*args) return Matrix([args]) elif isinstance(expr, Vector): - args = [[cls(_, evaluate=True)] for _ in expr.args[0]] + args = [[cls(a, evaluate=True)] for a in expr.args[0]] return ImmutableDenseMatrix(args) elif isinstance(expr, (Matrix, ImmutableDenseMatrix)): From 0ffdf1696d2ec60d0bfa4c529dc45b56d1ce2eee Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 10:39:33 +0100 Subject: [PATCH 13/21] rm commented line --- sympde/expr/evaluation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 0d2ae4fa..381c027e 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -620,7 +620,6 @@ def eval(cls, expr, domain): return expr.to_sympy() elif isinstance(expr, SympdeVector): -# return ImmutableDenseMatrix([expr.args[0]]) return expr.to_sympy() elif isinstance(expr, ScalarFunction): From 4a7d8c46fd0b8b4aa80d286f0b18c11b64cd4cf9 Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 10:43:19 +0100 Subject: [PATCH 14/21] fix exception for codacy --- sympde/core/matrices.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index 2f6ca676..44b80bdc 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -337,7 +337,7 @@ def __getitem__(self, key): try: i,j = key return self.args[0][i][j] - except: + except (ValueError,IndexError): return MatrixElement(self, key) def _sympystr(self, printer): @@ -389,8 +389,7 @@ def __getitem__(self, key): try: i = key return self.args[0][i] - except: - # TODO ARA do we keep returning a MatrixElement? + except (ValueError,IndexError): return MatrixElement(self, key) def _sympystr(self, printer): From 20e8c3141e0db9d1b173f5449bf3841522dbca99 Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 10:44:49 +0100 Subject: [PATCH 15/21] we should only name a variable as _ if it is a local variable which is not used later --- sympde/core/matrices.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index 44b80bdc..b55265da 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -406,8 +406,8 @@ def __eq__(self, a): return False def to_sympy(self): - _args = [[_] for _ in self.args[0]] - return ImmutableDenseMatrix(_args) + args = [[a] for a in self.args[0]] + return ImmutableDenseMatrix(args) Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = { From fb0c1e68f0ba7d9a2cb5bacd1684147ba106735e Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 10:50:43 +0100 Subject: [PATCH 16/21] rm try/except from __getitem__ --- sympde/core/matrices.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index b55265da..0e62a0b6 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -334,11 +334,8 @@ def name(self): return self._name def __getitem__(self, key): - try: - i,j = key - return self.args[0][i][j] - except (ValueError,IndexError): - return MatrixElement(self, key) + i,j = key + return self.args[0][i][j] def _sympystr(self, printer): sstr = printer.doprint @@ -386,11 +383,8 @@ def name(self): return self._name def __getitem__(self, key): - try: - i = key - return self.args[0][i] - except (ValueError,IndexError): - return MatrixElement(self, key) + i = key + return self.args[0][i] def _sympystr(self, printer): sstr = printer.doprint From 46b97b8a663f192c9390c31081bb6376b7399f83 Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 10:51:58 +0100 Subject: [PATCH 17/21] clean declaration of test functions --- sympde/expr/tests/test_expr.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 7dca172f..38cf049c 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -1914,8 +1914,8 @@ def test_matrices_vectors(): V = ScalarFunctionSpace('V', domain) W = VectorFunctionSpace('W', domain) - u,v = [element_of(V, name=i) for i in ['u', 'v']] - F,G = [element_of(W, name=i) for i in ['F', 'G']] + u, v = elements_of(V, names='u, v') + F, G = elements_of(W, names='F, G') #A = SympdeMatrix([[cos(theta), 0], [0, sin(theta)]], name='A') A = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A') From 93815b600c99ad42c363454a00516851d6d58488 Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 12:17:37 +0100 Subject: [PATCH 18/21] improve test equality between two Matrices or Vectors + add tests --- sympde/core/matrices.py | 36 ++++++++++++++++++++++++---------- sympde/expr/tests/test_expr.py | 28 ++++++++++++++++++++------ 2 files changed, 48 insertions(+), 16 deletions(-) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index 0e62a0b6..adf95d77 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -7,6 +7,8 @@ from sympy import ImmutableDenseMatrix +import functools + class MatrixSymbolicExpr(Expr): is_commutative = False @@ -344,11 +346,15 @@ def _sympystr(self, printer): def __hash__(self): return hash((self.name, self.args)) - def __eq__(self, a): - if isinstance(a, Matrix): - eq = self.name == a.name - return eq - return False + def __eq__(self, other): + if not isinstance(other, Matrix): + return False + + if not self.name == other.name: + return False + + result = all(x == y for x, y in zip(self.args[0], other.args[0])) + return result def to_sympy(self): return ImmutableDenseMatrix(self.args[0]) @@ -393,11 +399,21 @@ def _sympystr(self, printer): def __hash__(self): return hash((self.name, self.args)) - def __eq__(self, a): - if isinstance(a, Matrix): - eq = self.name == a.name - return eq - return False + def __eq__(self, other): + # TODO BUG + # We should check that other is a Vector + # right now, there is a problem when using Cross + # see the linearity test of + # f = lambda u,v: cross(b*u, b*v) + # where b is a Vector + if not isinstance(other, Vector): + return False + + if not self.name == other.name: + return False + + result = self.args[0] == other.args[0] + return result def to_sympy(self): args = [[a] for a in self.args[0]] diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 38cf049c..bed4464d 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -1921,6 +1921,22 @@ def test_matrices_vectors(): A = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A') b = Vector([1, 2, 3], name='b') + # ... + b1 = b + assert b1 == b + + b1 = Vector([1, 2, 3], name='b1') + assert not b1 == b + # ... + + # ... + A1 = A + assert A1 == A + + A1 = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A1') + assert not A1 == A + # ... + ######################################################################## # SCALAR-FUNCTION CASE ######################################################################## @@ -2005,12 +2021,12 @@ def test_matrices_vectors(): assert is_linear_expression(f(u,v), (v,)) # ... - # ... - f = lambda u,v: cross(b*u, b*v) - - assert is_linear_expression(f(u,v), (u,)) - assert is_linear_expression(f(u,v), (v,)) - # ... +# # ... +# f = lambda u,v: cross(b*u, b*v) +# +# assert is_linear_expression(f(u,v), (u,)) +# assert is_linear_expression(f(u,v), (v,)) +# # ... # ... f = lambda u,v: cross(b*u, grad(v)) From 7f83e33083edbae50e2b7b1663fd83f8d5b817aa Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 12:44:35 +0100 Subject: [PATCH 19/21] improve Vector & Matrix class constructors --- sympde/core/matrices.py | 73 ++++++++++++++--------------------------- 1 file changed, 24 insertions(+), 49 deletions(-) diff --git a/sympde/core/matrices.py b/sympde/core/matrices.py index adf95d77..4cfe7c30 100644 --- a/sympde/core/matrices.py +++ b/sympde/core/matrices.py @@ -296,44 +296,29 @@ def _sympystr(self, printer): return '{}[{}]'.format(sstr(self.args[0]),sstr(self.args[1])) class Matrix(MatrixSymbolicExpr): + def __new__(cls, mat, *, name): + if not isinstance(mat, list): + raise TypeError('Positional argument `mat` should be a list of lists.') - def __new__(cls, *args, **options): - name = options.pop('name') - if name is None: - raise ValueError('Expecting a name keyword.') - - if not( len(args) == 1 ): - raise ValueError('Expecting one argument.') - - _args = args[0] - if not isinstance(_args, list): - raise TypeError('Expecting a list.') - - for _ in _args: - if not isinstance(_, list): - raise TypeError('args components must be a list.') - - sizes = [len(_) for _ in _args] - if not( len(set(sizes)) == 1): - raise ValueError('Wrong matrix size') - - # make _args a Tuple of Tuples - newargs = [] - for _ in _args: - a = Tuple(*_) - newargs.append(a) - _args = Tuple(*newargs) + for row in mat: + if not isinstance(row, list): + raise TypeError('Each row of `mat` should be a list.') - args = (_args, args[1:]) + row_sizes = {len(row) for row in mat} + if len(row_sizes) != 1: + raise ValueError('Each row of `mat` should have the same length.') - obj = Expr.__new__(cls, *args) - obj._name = name + if not isinstance(name, str): + raise TypeError('Keyword-only argument `name` must be a string.') + # Call constructor of base class Expr with immutable arguments + new_mat = Tuple(*[Tuple(*row) for row in mat]) + obj = Expr.__new__(cls, new_mat, name) return obj @property def name(self): - return self._name + return self.args[1] def __getitem__(self, key): i,j = key @@ -362,31 +347,21 @@ def to_sympy(self): class Vector(MatrixSymbolicExpr): - def __new__(cls, *args, **options): - name = options.pop('name') - if name is None: - raise ValueError('Expecting a name keyword.') - - if not( len(args) == 1 ): - raise ValueError('Expecting one argument.') - - _args = args[0] - if not isinstance(_args, list): - raise TypeError('Expecting a list.') - - # make _args a Tuple - _args = Tuple(*_args) - - args = (_args, args[1:]) + def __new__(cls, vec, *, name): + if not isinstance(vec, list): + raise TypeError('Positional argument `vec` should be a list.') - obj = Expr.__new__(cls, *args) - obj._name = name + if not isinstance(name, str): + raise TypeError('Keyword-only argument `name` must be a string.') + # Call constructor of base class Expr with immutable arguments + new_vec = Tuple(*vec) + obj = Expr.__new__(cls, new_vec, name) return obj @property def name(self): - return self._name + return self.args[1] def __getitem__(self, key): i = key From 74e7475786e1fb5fde2968f6a9e7344ef7155afd Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 2 Mar 2024 13:10:13 +0100 Subject: [PATCH 20/21] Check equality of Cross arguments after simplifications --- sympde/calculus/core.py | 8 ++++++++ sympde/expr/tests/test_expr.py | 18 ++++++++++++------ 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index ac59ef41..151c0c01 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -428,6 +428,14 @@ def __new__(cls, arg1, arg2, **options): a,b = b,a c = -c + # after simplifications, the expression can be reduced to 0 in some cases + # here's an example; + # if b is a Vector and u and v are scalar functions, + # we should have the following + # cross(b*u,b*v) = u*v*cross(b,b) = 0 + if a == b: + return S.Zero + obj = Basic.__new__(cls, a, b) return c*obj diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index bed4464d..ae87fb58 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -1937,6 +1937,12 @@ def test_matrices_vectors(): assert not A1 == A # ... + # ... some tests on cross operator & vectors + assert cross(b,b) == 0 + assert cross(b*u,b*u) == 0 + assert cross(b*u,b*v) == 0 + # ... + ######################################################################## # SCALAR-FUNCTION CASE ######################################################################## @@ -2021,12 +2027,12 @@ def test_matrices_vectors(): assert is_linear_expression(f(u,v), (v,)) # ... -# # ... -# f = lambda u,v: cross(b*u, b*v) -# -# assert is_linear_expression(f(u,v), (u,)) -# assert is_linear_expression(f(u,v), (v,)) -# # ... + # ... + f = lambda u,v: cross(b*u, b*v) + + assert is_linear_expression(f(u,v), (u,)) + assert is_linear_expression(f(u,v), (v,)) + # ... # ... f = lambda u,v: cross(b*u, grad(v)) From 67d4e00fba2145bfd0f01ccdbf01d5d9e58018e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Sat, 2 Mar 2024 23:28:19 +0100 Subject: [PATCH 21/21] Clean some white spaces in test_expr.py --- sympde/expr/tests/test_expr.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index ae87fb58..0231b6a4 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -1882,8 +1882,8 @@ def test_interface_integral_4(): assert len(expr) == 1 assert expr[0].target == B.get_boundary(axis=0, ext=-1) - assert expr[0].expr[0,0] == u2[0]*v2[0] - assert expr[0].expr[1,1] == u2[1]*v2[1] + assert expr[0].expr[0,0] == u2[0]*v2[0] + assert expr[0].expr[1,1] == u2[1]*v2[1] a2 = BilinearForm((u2,v2), integral(I, dot(minus(u2),minus(v2)))) @@ -1891,12 +1891,10 @@ def test_interface_integral_4(): assert len(expr) == 1 assert expr[0].target == A.get_boundary(axis=0, ext=1) - assert expr[0].expr[0,0] == u2[0]*v2[0] - assert expr[0].expr[1,1] == u2[1]*v2[1] - + assert expr[0].expr[0,0] == u2[0]*v2[0] + assert expr[0].expr[1,1] == u2[1]*v2[1] # ... - #============================================================================== def test_matrices_vectors(): @@ -1905,7 +1903,7 @@ def test_matrices_vectors(): kappa = Constant('kappa', is_real=True) mu = Constant('mu', is_real=True) - theta = Constant('theta', is_real=True) + theta = Constant('theta', is_real=True) a1 = Constant('a1', is_real=True) a2 = Constant('a2', is_real=True) b1 = Constant('b1', is_real=True) @@ -1914,8 +1912,8 @@ def test_matrices_vectors(): V = ScalarFunctionSpace('V', domain) W = VectorFunctionSpace('W', domain) - u, v = elements_of(V, names='u, v') - F, G = elements_of(W, names='F, G') + u, v = elements_of(V, names='u, v') + F, G = elements_of(W, names='F, G') #A = SympdeMatrix([[cos(theta), 0], [0, sin(theta)]], name='A') A = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A') @@ -2110,7 +2108,6 @@ def test_matrices_vectors(): assert is_linear_expression(f(F,G), (G,)) # ... - #============================================================================== # CLEAN UP SYMPY NAMESPACE #==============================================================================