From 6d702d7061a0639c024c2dc1f12643b314a7c3c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 6 May 2025 12:23:39 +0200 Subject: [PATCH 01/15] Move 'eval' classmethod to InnerBasic from its subclasses --- sympde/core/algebra.py | 35 ++++++++--------------------------- 1 file changed, 8 insertions(+), 27 deletions(-) diff --git a/sympde/core/algebra.py b/sympde/core/algebra.py index 7b626117..476dc370 100644 --- a/sympde/core/algebra.py +++ b/sympde/core/algebra.py @@ -274,48 +274,29 @@ def __new__(cls, *args, **options): else: return r -class Inner_2d(InnerBasic): - @classmethod def eval(cls, *_args): - """.""" if not _args: return - if not( len(_args) == 2): + if len(_args) != 2: raise ValueError('Expecting two arguments') - u = _args[0] - v = _args[1] + u, v = _args u = Matrix(u) v = Matrix(v) # TODO add conjugate - M = u.transpose()*v + M = u.transpose() * v return M.trace() -class Inner_3d(InnerBasic): - - @classmethod - def eval(cls, *_args): - """.""" - - if not _args: - return - - if not( len(_args) == 2): - raise ValueError('Expecting two arguments') - - u = _args[0] - v = _args[1] - - u = Matrix(u) - v = Matrix(v) +# ... +class Inner_2d(InnerBasic): + pass - # TODO add conjugate - M = u.transpose()*v - return M.trace() # ... +class Inner_3d(InnerBasic): + pass From dbbea1b10064bc424c98303c047403e3f1fbb08c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 6 May 2025 12:24:34 +0200 Subject: [PATCH 02/15] Add Inner_1d class --- sympde/core/algebra.py | 3 +++ sympde/expr/evaluation.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/sympde/core/algebra.py b/sympde/core/algebra.py index 476dc370..bd2513d1 100644 --- a/sympde/core/algebra.py +++ b/sympde/core/algebra.py @@ -292,6 +292,9 @@ def eval(cls, *_args): M = u.transpose() * v return M.trace() +# ... +class Inner_1d(InnerBasic): + pass # ... class Inner_2d(InnerBasic): diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 9edda70a..85dbb039 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -16,7 +16,7 @@ from sympde.core.basic import _coeffs_registery from sympde.core.basic import CalculusFunction -from sympde.core.algebra import (Dot_1d, +from sympde.core.algebra import (Dot_1d, Inner_1d, Dot_2d, Inner_2d, Cross_2d, Dot_3d, Inner_3d, Cross_3d) from sympde.core.utils import random_string From 4e596397946e6e6c9625347354d0ef84f3f4777d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 11:04:55 +0200 Subject: [PATCH 03/15] Make ScalarFunction objects not commutative This allows writing inner(u,v) where (u,v) are scalar functions. --- sympde/topology/space.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/topology/space.py b/sympde/topology/space.py index cce1f556..45490d35 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -431,7 +431,7 @@ class ScalarFunction(Symbol): >>> V = SplineFemSpace('V') >>> phi = ScalarFunction(V, 'phi') """ - is_commutative = True + is_commutative = False _space = None _projection_of = None From 84d40f39699f73d49ee87bf717ca3d8e44a323f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 11:12:12 +0200 Subject: [PATCH 04/15] Convert scalar arguments of InnerBasic to SymPy Matrix This makes `TerminalExpr(inner(u, v), domain)` work, returning `u * v`. --- sympde/core/algebra.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/sympde/core/algebra.py b/sympde/core/algebra.py index bd2513d1..ceb1eaf9 100644 --- a/sympde/core/algebra.py +++ b/sympde/core/algebra.py @@ -285,6 +285,11 @@ def eval(cls, *_args): u, v = _args + # Matrix constructor requires a sequence + if not is_sequence(u): u = [u] + if not is_sequence(v): v = [v] + + # Convert arguments to Matrix (nothing happens if already Matrix) u = Matrix(u) v = Matrix(v) From 7309ff0db569f9138f0a0cd265c417ee5c1b40c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 12:14:47 +0200 Subject: [PATCH 05/15] Replace dot with inner in test_equation_2d.py --- sympde/expr/tests/test_equation_2d.py | 38 +++++++++++++-------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/sympde/expr/tests/test_equation_2d.py b/sympde/expr/tests/test_equation_2d.py index cb446f2f..95a3b425 100644 --- a/sympde/expr/tests/test_equation_2d.py +++ b/sympde/expr/tests/test_equation_2d.py @@ -50,7 +50,7 @@ def test_equation_2d_1(): int_2 = lambda expr: integral(B2, expr) # ... bilinear/linear forms - expr = dot(grad(v), grad(u)) + expr = inner(grad(v), grad(u)) a1 = BilinearForm((v,u), int_0(expr)) expr = v*u @@ -121,12 +121,12 @@ def test_equation_2d_1(): equation = Equation(a1, l1, tests=v, trials=u, bc=EssentialBC(u, 0, B1)) # ... - # ... Poisson with Nitsch method - g = cos(pi*x)*cos(pi*y) - a0 = BilinearForm((u,v), int_0(dot(grad(u),grad(v)))) - a_B1 = BilinearForm((u,v), - int_1(kappa * u*dot(grad(v), nn) - - v*dot(grad(u), nn) - + u*v/eps)) + # ... Poisson with Nitsche's method + g = cos(pi*x) * cos(pi*y) + a0 = BilinearForm((u, v), int_0(inner(grad(u), grad(v)))) + a_B1 = BilinearForm((u, v), - int_1(kappa * u * dot(grad(v), nn) + - v * dot(grad(u), nn) + + u * v/eps)) a = BilinearForm((u,v), a0(u,v) + a_B1(u,v)) @@ -170,12 +170,12 @@ def test_equation_2d_2(): int_0 = lambda expr: integral(domain , expr) - s = BilinearForm((tau,sigma), int_0(dot(grad(tau), grad(sigma)))) + s = BilinearForm((tau,sigma), int_0(inner(grad(tau), grad(sigma)))) m = BilinearForm((tau,sigma), int_0(tau*sigma)) b1 = BilinearForm((tau,dw), int_0(bracket(pn, dw) * tau)) b2 = BilinearForm((tau,dp), int_0(bracket(dp, wn) * tau)) - l1 = LinearForm(tau, int_0(bracket(pn, wn)*tau - 1./Re * dot(grad(tau), grad(wn)))) + l1 = LinearForm(tau, int_0(bracket(pn, wn)*tau - 1./Re * inner(grad(tau), grad(wn)))) expr = m(tau,dw) - alpha*dt*b1(tau,dw) - dt*b2(tau,dp) - (alpha*dt/Re)*s(tau,dw) a = BilinearForm(((tau, sigma),(dp,dw)), expr) @@ -207,7 +207,7 @@ def test_equation_2d_3(): int_1 = lambda expr: integral(B1, expr) # ... bilinear/linear forms - a1 = BilinearForm((v,u), int_0(dot(grad(v), grad(u)))) + a1 = BilinearForm((v,u), int_0(inner(grad(v), grad(u)))) a2 = BilinearForm((v,u), int_0(v*u)) l1 = LinearForm(v, int_0(x*y*v)) @@ -243,7 +243,7 @@ def test_equation_2d_4(): a1 = BilinearForm((v,u), int_0(inner(grad(v), grad(u)))) f = Matrix([x*y, sin(pi*x)*sin(pi*y)]) - l1 = LinearForm(v, int_0(dot(f,v))) + l1 = LinearForm(v, int_0(inner(f, v))) # ... # ... @@ -291,7 +291,7 @@ def test_equation_2d_5(): a = BilinearForm(((v,q),(u,p)), a0(v,u) + a1(q,p)) print(' a done.') - l0 = LinearForm(v, int_0(dot(f0, v))) + l0 = LinearForm(v, int_0(inner(f0, v))) l1 = LinearForm(q, int_0(f1*q)) l = LinearForm((v,q), l0(v) + l1(q)) @@ -337,7 +337,7 @@ def test_equation_2d_6(): int_0 = lambda expr: integral(domain , expr) # ... - expr = kappa * dot(grad(u), grad(v)) + dot(b, grad(u)) * v + expr = kappa * inner(grad(u), grad(v)) + inner(b, grad(u)) * v a = BilinearForm((v,u), int_0(expr)) # ... # ... @@ -346,26 +346,26 @@ def test_equation_2d_6(): # ... # ... - expr = (- kappa * laplace(u) + dot(b, grad(u))) * dot(b, grad(v)) + expr = (- kappa * laplace(u) + inner(b, grad(u))) * inner(b, grad(v)) s1 = BilinearForm((v,u), int_0(expr)) - expr = - f * dot(b, grad(v)) + expr = - f * inner(b, grad(v)) l1 = LinearForm(v, int_0(expr)) # ... # ... - expr = (- kappa * laplace(u) + dot(b, grad(u))) * ( dot(b, grad(v)) - kappa * laplace(v)) + expr = (- kappa * laplace(u) + inner(b, grad(u))) * ( inner(b, grad(v)) - kappa * laplace(v)) s2 = BilinearForm((v,u), int_0(expr)) - expr = - f * ( dot(b, grad(v)) - kappa * laplace(v)) + expr = - f * ( inner(b, grad(v)) - kappa * laplace(v)) l2 = LinearForm(v, int_0(expr)) # ... # ... - expr = (- kappa * laplace(u) + dot(b, grad(u))) * ( dot(b, grad(v)) + kappa * laplace(v)) + expr = (- kappa * laplace(u) + inner(b, grad(u))) * ( inner(b, grad(v)) + kappa * laplace(v)) s3 = BilinearForm((v,u), int_0(expr)) - expr = - f * ( dot(b, grad(v)) + kappa * laplace(v)) + expr = - f * ( inner(b, grad(v)) + kappa * laplace(v)) l3 = LinearForm(v, int_0(expr)) # ... From 4887e7cf2355058739b220292c19524b7b129fa6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 12:20:52 +0200 Subject: [PATCH 06/15] Revert "Make ScalarFunction objects not commutative" This reverts commit 4e596397946e6e6c9625347354d0ef84f3f4777d. --- sympde/topology/space.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/topology/space.py b/sympde/topology/space.py index 45490d35..cce1f556 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -431,7 +431,7 @@ class ScalarFunction(Symbol): >>> V = SplineFemSpace('V') >>> phi = ScalarFunction(V, 'phi') """ - is_commutative = False + is_commutative = True _space = None _projection_of = None From d5cba77a334afb3749a2d49ca23b93599fe11427 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 13:07:34 +0200 Subject: [PATCH 07/15] Test Inner constructor with ScalarFunction arguments --- sympde/calculus/tests/test_calculus.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sympde/calculus/tests/test_calculus.py b/sympde/calculus/tests/test_calculus.py index f7ea8fac..0af8ac55 100644 --- a/sympde/calculus/tests/test_calculus.py +++ b/sympde/calculus/tests/test_calculus.py @@ -76,11 +76,14 @@ def test_Cross(dim): # TODO: check exceptions #============================================================================== +@pytest.mark.parametrize('vector', [True, False]) @pytest.mark.parametrize('dim', [1, 2, 3]) -def test_Inner(dim): +def test_Inner(dim, vector): + + FunctionSpace = VectorFunctionSpace if vector else ScalarFunctionSpace domain = Domain('Omega', dim=dim) - W = VectorFunctionSpace('W', domain) + W = FunctionSpace('W', domain) a, b, c = elements_of(W, names='a, b, c') r = Constant('r') From efe736c68216f33d3ad374b790be6843320a8c55 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 13:11:24 +0200 Subject: [PATCH 08/15] Make Inner constructor accept ScalarFunction arguments --- sympde/calculus/core.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index ac59ef41..0eb21419 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -500,25 +500,27 @@ def __new__(cls, arg1, arg2, **options): else: b = [arg2] - args_1 = [i for i in a if not i.is_commutative] + args_1 = [i for i in a if not i.is_commutative or has(i, types)] c1 = [i for i in a if not i in args_1] - args_2 = [i for i in b if not i.is_commutative] + args_2 = [i for i in b if not i.is_commutative or has(i, types)] c2 = [i for i in b if not i in args_2] a = reduce(mul, args_1) b = reduce(mul, args_2) - c = Mul(*c1)*Mul(*c2) + c = Mul(*c1) * Mul(*c2) + # TODO [YG 20.05.2025] : apply conjugate if arguments are flipped if str(a) > str(b): - a,b = b,a + a, b = b, a obj = Basic.__new__(cls, a, b) + # This is a norm when a == b if a == b: obj.is_real = True obj.is_positive = True - return c*obj + return c * obj #============================================================================== class Outer(BasicOperator): From 3ca0f185ca53d442ece97e65508d7754149ebed0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 15:58:38 +0200 Subject: [PATCH 09/15] Check TerminalExpr output in test_expr.py::test_terminal_expr_linear_2d_1 --- sympde/expr/tests/test_expr.py | 94 ++++++++++++++++++++++------------ 1 file changed, 62 insertions(+), 32 deletions(-) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index c14fe9be..2205ff89 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -32,7 +32,7 @@ from sympde.expr.expr import integral from sympde.expr.expr import Functional, Norm from sympde.expr.expr import linearize -from sympde.expr.evaluation import TerminalExpr +from sympde.expr.evaluation import TerminalExpr, DomainExpression, BoundaryExpression #============================================================================== @@ -369,7 +369,8 @@ def test_terminal_expr_linear_2d_1(): domain = Domain('Omega', dim=2) B1 = Boundary(r'\Gamma_1', domain) - x,y = domain.coordinates + # Domain has no mapping, so physical coordinates are the same as logical + x, y = x1, x2 = domain.coordinates kappa = Constant('kappa', is_real=True) mu = Constant('mu' , is_real=True) @@ -380,69 +381,98 @@ def test_terminal_expr_linear_2d_1(): 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']] - # ... int_0 = lambda expr: integral(domain , expr) int_1 = lambda expr: integral(B1, expr) - l = LinearForm(v, int_0(x*y*v)) + # ... + expr = x * y * v + l = LinearForm(v, int_0(expr)) + l_term = TerminalExpr(l, domain) - print(TerminalExpr(l, domain)) - print('') + l_term_expected = (DomainExpression(domain.interior, expr),) + assert l_term == l_term_expected # ... # ... - l = LinearForm(v, int_0(x*y*v + v)) - print(TerminalExpr(l, domain)) - print('') + expr = x * y * v + v + l = LinearForm(v, int_0(expr)) + l_term = TerminalExpr(l, domain) + + l_term_expected = (DomainExpression(domain.interior, expr),) + assert l_term == l_term_expected # ... # ... g = Matrix((x**2, y**2)) - l = LinearForm(v, int_1(v*dot(g, nn))) - print(TerminalExpr(l, domain)) - print('') + expr = v * dot(g, nn) + l = LinearForm(v, int_1(expr)) + l_term = TerminalExpr(l, domain) + + expr_2 = v * (g[0] * nn[0] + g[1] * nn[1]) + l_term_expected = (BoundaryExpression(B1, expr_2),) + assert l_term == l_term_expected # ... # ... g = Matrix((x**2, y**2)) - l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v)) - print(TerminalExpr(l, domain)) - print('') + expr_0 = x * y * v + expr_1 = v * dot(g, nn) + l = LinearForm(v, int_1(expr_1) + int_0(expr_0)) + l_term = TerminalExpr(l, domain) + + expr_2 = v * (g[0] * nn[0] + g[1] * nn[1]) + l_term_expected = (DomainExpression(domain.interior, expr_0), + BoundaryExpression(B1, expr_2)) + assert l_term == l_term_expected # ... # ... - l1 = LinearForm(v1, int_0(x*y*v1)) + l1 = LinearForm(v1, int_0(x * y * v1)) l = LinearForm(v, l1(v)) - print(TerminalExpr(l, domain)) - print('') + l_term = TerminalExpr(l, domain) + + l_term_expected = (DomainExpression(domain.interior, x * y * v),) + assert l_term == l_term_expected # ... # ... - g = Matrix((x,y)) - l1 = LinearForm(v1, int_0(x*y*v1)) + g = Matrix((x, y)) + l1 = LinearForm(v1, int_0(x * y * v1)) l2 = LinearForm(v2, int_0(dot(grad(v2), g))) + l = LinearForm(v, l1(v) + l2(v)) + l_term = TerminalExpr(l, domain) - l = LinearForm(v, l1(v) + l2(v)) - print(TerminalExpr(l, domain)) - print('') + expr_2 = v * x * y + x * dx1(v) + y * dx2(v) + l_term_expected = (DomainExpression(domain.interior, expr_2),) + assert l_term == l_term_expected # ... + # ... - l1 = LinearForm(v1, int_0(x*y*v1)) + l1 = LinearForm(v1, int_0(x * y * v1)) l2 = LinearForm(v1, int_0(v1)) - l = LinearForm(v, l1(v) + kappa*l2(v)) - print(TerminalExpr(l, domain)) - print('') + l = LinearForm(v, l1(v) + kappa * l2(v)) + + l_term = TerminalExpr(l, domain) + + expr_2 = kappa * v + v * x * y + l_term_expected = (DomainExpression(domain.interior, expr_2),) + assert l_term == l_term_expected # ... # ... g = Matrix((x**2, y**2)) - l1 = LinearForm(v1, int_0(x*y*v1)) + l1 = LinearForm(v1, int_0(x * y * v1)) l2 = LinearForm(v1, int_0(v1)) - l3 = LinearForm(v, int_1(v*dot(g, nn))) - l = LinearForm(v, l1(v) + kappa*l2(v) + mu*l3(v)) - print(TerminalExpr(l, domain)) - print('') + l3 = LinearForm(v, int_1(v * dot(g, nn))) + l = LinearForm(v, l1(v) + kappa * l2(v) + mu * l3(v)) + l_term = TerminalExpr(l, domain) + + expr_2 = kappa * v + v * x * y + expr_3 = mu * v * (g[0] * nn[0] + g[1] * nn[1]) + l_term_expected = (DomainExpression(domain.interior, expr_2), + BoundaryExpression(B1, expr_3)) + assert l_term == l_term_expected # ... #============================================================================== From 797704eb911af7b8ec536e55c9af5fc0d539f3de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 17:22:05 +0200 Subject: [PATCH 10/15] Remove unused imports from core/algebra.py --- sympde/core/algebra.py | 23 +---------------------- 1 file changed, 1 insertion(+), 22 deletions(-) diff --git a/sympde/core/algebra.py b/sympde/core/algebra.py index ceb1eaf9..bf83ec04 100644 --- a/sympde/core/algebra.py +++ b/sympde/core/algebra.py @@ -1,38 +1,17 @@ # coding: utf-8 -import numpy as np -from itertools import groupby - #from sympy.core.sympify import sympify from sympy.simplify.simplify import simplify -from sympy import Symbol -from sympy import Lambda -from sympy import Function -from sympy import bspline_basis -from sympy import lambdify -from sympy import cos -from sympy import sin -from sympy import Rational -from sympy import diff from sympy import Matrix, ImmutableDenseMatrix -from sympy import latex -from sympy import I as sympy_I from sympy.core import Basic from sympy.core.singleton import S -from sympy.simplify.simplify import nsimplify -from sympy.utilities.lambdify import implemented_function -from sympy.matrices.dense import MutableDenseMatrix from sympy import Mul, Add -from sympy import postorder_traversal -from sympy import preorder_traversal -from sympy.core.expr import Expr from sympy.core.containers import Tuple -from sympy import Integer, Float from sympy import Add, Mul from sympy import simplify from sympy import S from sympy import Basic -from sympy import Indexed, IndexedBase +from sympy import Indexed from sympde.old_sympy_utilities import is_sequence from .basic import CalculusFunction From e920d3d3b05802fac687109eddc3571e965c8426 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 19:50:49 +0200 Subject: [PATCH 11/15] TerminalExpr: replace string evaluations w/ class dictionaries --- sympde/expr/evaluation.py | 76 +++++++++++++++++++++++++++++++-------- 1 file changed, 61 insertions(+), 15 deletions(-) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 85dbb039..426058e0 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -49,15 +49,22 @@ from sympde.topology.derivatives import get_atom_logical_derivatives from sympde.topology.derivatives import dx, dy, dz from sympde.topology.derivatives import dx1, dx2, dx3 -from sympde.topology.derivatives import (Grad_1d, Div_1d, - Grad_2d, Curl_2d, Rot_2d, Div_2d, - Grad_3d, Curl_3d, Div_3d) + +from sympde.calculus.core import Grad, Div, Curl, Rot, Bracket, Laplace, Hessian +from sympde.calculus.core import Dot, Inner, Cross + +from sympde.topology.derivatives import Grad_1d, Grad_2d, Grad_3d +from sympde.topology.derivatives import Div_1d, Div_2d, Div_3d +from sympde.topology.derivatives import Curl_2d, Curl_3d +from sympde.topology.derivatives import Rot_2d from sympde.topology.derivatives import Bracket_2d from sympde.topology.derivatives import Laplace_1d, Laplace_2d, Laplace_3d from sympde.topology.derivatives import Hessian_1d, Hessian_2d, Hessian_3d -from sympde.topology.derivatives import (LogicalGrad_1d, LogicalDiv_1d, - LogicalGrad_2d, LogicalCurl_2d, LogicalRot_2d, LogicalDiv_2d, - LogicalGrad_3d, LogicalCurl_3d, LogicalDiv_3d) + +from sympde.topology.derivatives import LogicalGrad_1d, LogicalGrad_2d, LogicalGrad_3d +from sympde.topology.derivatives import LogicalDiv_1d, LogicalDiv_2d, LogicalDiv_3d +from sympde.topology.derivatives import LogicalCurl_2d, LogicalCurl_3d +from sympde.topology.derivatives import LogicalRot_2d from sympde.topology.derivatives import LogicalBracket_2d from sympde.topology.derivatives import LogicalLaplace_1d, LogicalLaplace_2d, LogicalLaplace_3d from sympde.topology.derivatives import LogicalHessian_1d, LogicalHessian_2d, LogicalHessian_3d @@ -91,6 +98,33 @@ 'is_sequence', ) +#============================================================================== +differential_operators = { + Grad : (None, Grad_1d, Grad_2d, Grad_3d), + Div : (None, Div_1d, Div_2d, Div_3d), + Curl : (None, None, Curl_2d, Curl_3d), + Rot : (None, None, Rot_2d, None), + Bracket: (None, None, Bracket_2d, None), + Laplace: (None, Laplace_1d, Laplace_2d, Laplace_3d), + Hessian: (None, Hessian_1d, Hessian_2d, Hessian_3d), +} + +logical_differential_operators = { + Grad : (None, LogicalGrad_1d, LogicalGrad_2d, LogicalGrad_3d), + Div : (None, LogicalDiv_1d, LogicalDiv_2d, LogicalDiv_3d), + Curl : (None, None, LogicalCurl_2d, LogicalCurl_3d), + Rot : (None, None, LogicalRot_2d, None), + Bracket: (None, None, LogicalBracket_2d, None), + Laplace: (None, LogicalLaplace_1d, LogicalLaplace_2d, LogicalLaplace_3d), + Hessian: (None, LogicalHessian_1d, LogicalHessian_2d, LogicalHessian_3d), +} + +generic_operators = { + Dot : (None, Dot_1d, Dot_2d, Dot_3d), + Inner: (None, Inner_1d, Inner_2d, Inner_3d), + Cross: (None, None, Cross_2d, Cross_3d), +} + #============================================================================== def is_sequence(a): return isinstance(a, (list,tuple,Tuple)) @@ -769,23 +803,35 @@ def eval(cls, expr, domain): elif isinstance(expr, BasicExpr): return cls.eval(expr.expr, domain=domain) - elif isinstance(expr, _diff_ops): - op = type(expr) + elif isinstance(expr, (*differential_operators,)): if domain.mapping is None: - new = eval('Logical{0}_{1}d'.format(op, dim)) + new = logical_differential_operators[type(expr)][dim] else: - new = eval('{0}_{1}d'.format(op, dim)) - + new = differential_operators[type(expr)][dim] args = [cls.eval(i, domain=domain) for i in expr.args] return new(*args) - elif isinstance(expr, _generic_ops): - # if i = Dot(...) then type(i) is Grad - op = type(expr) - new = eval('{0}_{1}d'.format(op, dim)) + elif isinstance(expr, (*generic_operators,)): + new = generic_operators[type(expr)][dim] args = [cls.eval(i, domain=domain) for i in expr.args] return new(*args) +# elif isinstance(expr, _diff_ops): +# op = type(expr) +# if domain.mapping is None: +# new = eval('Logical{0}_{1}d'.format(op, dim)) +# else: +# new = eval('{0}_{1}d'.format(op, dim)) +# args = [cls.eval(i, domain=domain) for i in expr.args] +# return new(*args) +# +# elif isinstance(expr, _generic_ops): +# # if i = Dot(...) then type(i) is Grad +# op = type(expr) +# new = eval('{0}_{1}d'.format(op, dim)) +# args = [cls.eval(i, domain=domain) for i in expr.args] +# return new(*args) + elif isinstance(expr, Trace): # TODO treate different spaces if expr.order == 0: From 593a4a3a42a14e6c68e019df05aad492862721e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 20 May 2025 19:52:59 +0200 Subject: [PATCH 12/15] Remove unused imports from expr/evaluation.py --- sympde/expr/evaluation.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 426058e0..40d57718 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -1,13 +1,7 @@ # coding: utf-8 -from itertools import product - - -import numpy as np - from sympy import Abs, S, cacheit from sympy import Indexed, Matrix, ImmutableDenseMatrix -from sympy import expand from sympy.core import Basic, Symbol from sympy.core import Add, Mul, Pow from sympy.core.expr import AtomicExpr @@ -19,19 +13,17 @@ from sympde.core.algebra import (Dot_1d, Inner_1d, Dot_2d, Inner_2d, Cross_2d, Dot_3d, Inner_3d, Cross_3d) -from sympde.core.utils import random_string -from sympde.calculus import jump, avg, minus, plus +from sympde.calculus import jump, 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.calculus.matrices import MatrixElement, SymbolicTrace from sympde.topology.basic import BasicDomain, Union, Interval from sympde.topology.basic import Boundary, Interface from sympde.topology.basic import InteriorDomain from sympde.topology.domain import NormalVector, TangentVector, NCube, NCubeInterior -from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, MultiPatchMapping, JacobianInverseSymbol +from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, JacobianInverseSymbol from sympde.topology.mapping import LogicalExpr, PullBack # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping @@ -40,7 +32,6 @@ from sympde.topology.space import VectorFunction from sympde.topology.space import IndexedVectorFunction from sympde.topology.space import Trace -from sympde.topology.space import element_of from sympde.topology.space import ScalarFunctionSpace from sympde.topology.derivatives import _partial_derivatives From 923e2c9e21c6df3efa6d495b77a52bcd99b41089 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 22 May 2025 20:07:07 +0200 Subject: [PATCH 13/15] Check TerminalExpr output in test_expr.py::test_terminal_expr_linear_2d_2 --- sympde/expr/tests/test_expr.py | 184 +++++++++++++++++++++++---------- 1 file changed, 128 insertions(+), 56 deletions(-) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 2205ff89..b194c1c5 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -481,10 +481,11 @@ def test_terminal_expr_linear_2d_2(): domain = Domain('Omega', dim=2) B1 = Boundary(r'\Gamma_1', domain) - x,y = domain.coordinates + x1, x2 = domain.coordinates kappa = Constant('kappa', is_real=True) mu = Constant('mu' , is_real=True) + nn = NormalVector('nn') V = VectorFunctionSpace('V', domain) @@ -495,68 +496,139 @@ 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 = Matrix((x1, x2)) l = LinearForm(v, int_0(dot(g, v))) - print(TerminalExpr(l, domain)) - print('') + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([[x1 * v[0]], [x2 * v[1]]]) + l_term_expected = (DomainExpression(domain.interior, expr_0),) + assert l_term == l_term_expected # ... # ... - g = Matrix((x,y)) + g = Matrix((x1, x2)) l = LinearForm(v, int_0(dot(g, v) + div(v))) - print(TerminalExpr(l, domain)) - print('') + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([[x1 * v[0] + dx1(v[0])], [x2 * v[1] + dx2(v[1])]]) + l_term_expected = (DomainExpression(domain.interior, expr_0),) + assert l_term == l_term_expected # ... - # TODO -# # ... -# g = Tuple(x**2, y**2) -# l = LinearForm(v, v*trace_1(g, B1)) -# print(TerminalExpr(l)) -# print('') -# # ... -# -# # ... -# g = Tuple(x**2, y**2) -# l = LinearForm(v, v*trace_1(g, B1) + x*y*v) -# print(TerminalExpr(l)) -# print('') -# # ... -# -# # ... -# l1 = LinearForm(v1, x*y*v1) -# l = LinearForm(v, l1(v)) -# print(TerminalExpr(l)) -# print('') -# # ... -# -# # ... -# g = Tuple(x,y) -# l1 = LinearForm(v1, x*y*v1) -# l2 = LinearForm(v2, dot(grad(v2), g)) -# -# l = LinearForm(v, l1(v) + l2(v)) -# print(TerminalExpr(l)) -# print('') -# # ... -# -# # ... -# l1 = LinearForm(v1, x*y*v1) -# l2 = LinearForm(v1, v1) -# l = LinearForm(v, l1(v) + kappa*l2(v)) -# print(TerminalExpr(l)) -# print('') -# # ... -# -# # ... -# g = Tuple(x**2, y**2) -# l1 = LinearForm(v1, x*y*v1) -# l2 = LinearForm(v1, v1) -# l3 = LinearForm(v, v*trace_1(g, B1)) -# l = LinearForm(v, l1(v) + kappa*l2(v) + mu*l3(v)) -# print(TerminalExpr(l)) -# print('') -# # ... + # ... + g = Tuple(x1**2, x2**2) + l = LinearForm(v, int_1(v * dot(g, nn))) + l_term = TerminalExpr(l, domain) + + expr_1 = Matrix([ + [(x1**2 * nn[0] + x2**2 * nn[1]) * v[0]], + [ 0], + [ 0], + [(x1**2 * nn[0] + x2**2 * nn[1]) * v[1]] + ]) + l_term_expected = (BoundaryExpression(B1, expr_1),) + assert l_term == l_term_expected + # ... + + # ... + g = Tuple(x1**2, x2**2) + l = LinearForm(v, int_1(v * dot(g, nn)) + int_0(x1 * x2 * v)) + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([ + [x1 * x2 * v[0]], + [ 0], + [ 0], + [x1 * x2 * v[1]] + ]) + expr_1 = Matrix([ + [(x1**2 * nn[0] + x2**2 * nn[1]) * v[0]], + [ 0], + [ 0], + [(x1**2 * nn[0] + x2**2 * nn[1]) * v[1]] + ]) + l_term_expected = (DomainExpression(domain.interior, expr_0), + BoundaryExpression(B1, expr_1)) + assert l_term == l_term_expected + # ... + + # ... + l1 = LinearForm(v1, int_0(x1 * x2 * v1)) + l = LinearForm(v, l1(v)) + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([ + [x1 * x2 * v[0]], + [ 0], + [ 0], + [x1 * x2 * v[1]] + ]) + l_term_expected = (DomainExpression(domain.interior, expr_0),) + assert l_term == l_term_expected + # ... + + # ... + g = Tuple(x1, x2) + l1 = LinearForm(v1, int_0(x1 * x2 * v1)) + l2 = LinearForm(v2, int_1(dot(grad(v2), g))) + l = LinearForm(v, l1(v) + l2(v)) + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([ + [x1 * x2 * v[0]], + [ 0], + [ 0], + [x1 * x2 * v[1]] + ]) + expr_1 = Matrix([ + [x1 * dx1(v[0])], + [x2 * dx1(v[1])] + ]) + l_term_expected = (DomainExpression(domain.interior, expr_0), + BoundaryExpression(B1, expr_1)) + assert l_term == l_term_expected + # ... + + # ... + l1 = LinearForm(v1, int_0(x1 * x2 * v1)) + l2 = LinearForm(v1, int_0(v1)) + l = LinearForm(v, l1(v) + kappa * l2(v)) + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([ + [kappa * v[0] + x1 * x2 * v[0]], + [ 0], + [ 0], + [kappa * v[1] + x1 * x2 * v[1]] + ]) + l_term_expected = (DomainExpression(domain.interior, expr_0),) + assert l_term == l_term_expected + # ... + + # ... + g = Tuple(x1**2, x2**2) + l1 = LinearForm(v1, int_0(x1 * x2 * v1)) + l2 = LinearForm(v1, int_0(v1)) + l3 = LinearForm(v, int_1(v * dot(g, nn))) + l = LinearForm(v, l1(v) + kappa * l2(v) + mu * l3(v)) + l_term = TerminalExpr(l, domain) + + expr_0 = Matrix([ + [kappa * v[0] + x1 * x2 * v[0]], + [ 0], + [ 0], + [kappa * v[1] + x1 * x2 * v[1]] + ]) + expr_1 = Matrix([ + [mu * (x1**2 * nn[0] + x2**2 * nn[1]) * v[0]], + [ 0], + [ 0], + [mu * (x1**2 * nn[0] + x2**2 * nn[1]) * v[1]] + ]) + l_term_expected = (DomainExpression(domain.interior, expr_0), + BoundaryExpression(B1, expr_1)) + assert l_term == l_term_expected + # ... #============================================================================== def test_terminal_expr_linear_2d_3(): From 6a1841579c5d0a2f136cffa32ebbfd70366539e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 5 Jun 2025 19:33:17 +0200 Subject: [PATCH 14/15] Use inner for scalar inputs in function test_equation_2d_1 --- sympde/expr/tests/test_equation_2d.py | 32 +++++++++++++-------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/sympde/expr/tests/test_equation_2d.py b/sympde/expr/tests/test_equation_2d.py index 95a3b425..6abf4675 100644 --- a/sympde/expr/tests/test_equation_2d.py +++ b/sympde/expr/tests/test_equation_2d.py @@ -51,21 +51,21 @@ def test_equation_2d_1(): # ... bilinear/linear forms expr = inner(grad(v), grad(u)) - a1 = BilinearForm((v,u), int_0(expr)) + a1 = BilinearForm((u, v), int_0(expr)) - expr = v*u - a2 = BilinearForm((v,u),int_0(expr)) + expr = inner(u, v) + a2 = BilinearForm((u, v),int_0(expr)) - expr = v*u - a_B1 = BilinearForm((v, u), int_1(expr)) + expr = inner(u, v) + a_B1 = BilinearForm((u, v), int_1(expr)) - expr = x*y*v + expr = inner(x * y, v) l1 = LinearForm(v, int_0(expr)) - expr = cos(x+y)*v + expr = inner(cos(x + y), v) l2 = LinearForm(v, int_0(expr)) - expr = x*y*v + expr = inner(x * y, v) l_B2 = LinearForm(v, int_2(expr)) # ... @@ -96,23 +96,23 @@ def test_equation_2d_1(): # ... # ... - a = BilinearForm((v,u), a1(v,u) + alpha*a2(v,u)) + a = BilinearForm((u, v), a1(u, v) + alpha*a2(u, v)) equation = Equation(a, l1, tests=v, trials=u) # ... # ... - a = BilinearForm((v,u), a1(v,u) + a_B1(v,u)) + a = BilinearForm((u, v), a1(u, v) + a_B1(u, v)) equation = Equation(a, l1, tests=v, trials=u) # ... # ... - a = BilinearForm((v,u), a1(v,u) + a_B1(v,u)) + a = BilinearForm((u, v), a1(u, v) + a_B1(u, v)) l = LinearForm(v, l1(v) + l2(v)) equation = Equation(a, l, tests=v, trials=u) # ... # ... - a = BilinearForm((v,u), a1(v,u) + a_B1(v,u)) + a = BilinearForm((u, v), a1(u, v) + a_B1(u, v)) l = LinearForm(v, l1(v) + alpha*l_B2(v)) equation = Equation(a, l, tests=v, trials=u) # ... @@ -125,12 +125,12 @@ def test_equation_2d_1(): g = cos(pi*x) * cos(pi*y) a0 = BilinearForm((u, v), int_0(inner(grad(u), grad(v)))) a_B1 = BilinearForm((u, v), - int_1(kappa * u * dot(grad(v), nn) - - v * dot(grad(u), nn) - + u * v/eps)) + - inner(dot(grad(u), nn), v) + + inner(u, v)/eps)) - a = BilinearForm((u,v), a0(u,v) + a_B1(u,v)) - l = LinearForm(v, int_0(g*v/eps) - int_1(kappa*v*g)) + a = BilinearForm((u, v), a0(u, v) + a_B1(u, v)) + l = LinearForm(v, int_0(inner(g, v)) / eps - kappa * int_1(inner(g, v))) equation = Equation(a, l, tests=v, trials=u) # ... From a3848dd97052e099da4ad5214f82e51df621b6e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 5 Jun 2025 19:35:21 +0200 Subject: [PATCH 15/15] Fix bug in Inner.__new__, allowing inner w/ scalar inputs --- sympde/calculus/core.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index 0eb21419..593dd902 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -505,8 +505,8 @@ def __new__(cls, arg1, arg2, **options): args_2 = [i for i in b if not i.is_commutative or has(i, types)] c2 = [i for i in b if not i in args_2] - a = reduce(mul, args_1) - b = reduce(mul, args_2) + a = reduce(mul, args_1, S.One) + b = reduce(mul, args_2, S.One) c = Mul(*c1) * Mul(*c2) # TODO [YG 20.05.2025] : apply conjugate if arguments are flipped