diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index e75189be..9edda70a 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -567,7 +567,7 @@ def eval(cls, expr, domain): if isinstance(domain, (NCube, NCubeInterior)): bounds = domain.min_coords if ext == -1 else domain.max_coords J = J.subs(coordinates[axis], bounds[axis]) - newcoords = [Symbol(c.name+"_plus") for c in coordinates] + newcoords = [Symbol(c.name+"_plus", real=True) for c in coordinates] subs = list(zip(coordinates, newcoords)) J = J.subs(subs) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index ebcfb19b..c14fe9be 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -35,6 +35,7 @@ from sympde.expr.evaluation import TerminalExpr #============================================================================== + def test_linear_expr_2d_1(): domain = Domain('Omega', dim=2) @@ -897,6 +898,30 @@ def test_terminal_expr_bilinear_3d_1(): assert e2[0].expr == dx(um)*dx(vm) + dy(um)*dy(vm) + dz(um)*dz(vm) assert e3[0].expr.factor() == (dx1(u)*dx1(v) + dx2(u)*dx2(v) + dx3(u)*dx3(v))*det +#============================================================================== +def test_terminal_expr_1d_1(): + domain = Domain('Omega', dim=1) + x = domain.coordinates + assert TerminalExpr(grad(x), domain) == 1 + assert TerminalExpr(grad(exp(1j*x) - x + 7), domain) == 1j*exp(1j*x) - 1 + +#============================================================================== +def test_terminal_expr_2d_1(): + domain = Domain('Omega', dim=2) + x, y = domain.coordinates + assert TerminalExpr(grad(x), domain) == Matrix([[1], [0]]) + assert TerminalExpr(grad(y), domain) == Matrix([[0], [1]]) + assert TerminalExpr(grad(exp(1j*x)*sin(2*y)), domain) == Matrix([[1j*exp(1j*x)*sin(2*y)], [2*exp(1j*x)*cos(2*y)]]) + +#============================================================================== +def test_terminal_expr_3d_1(): + domain = Domain('Omega', dim=3) + x, y, z = domain.coordinates + assert TerminalExpr(grad(x), domain) == Matrix([[1], [0], [0]]) + assert TerminalExpr(grad(y), domain) == Matrix([[0], [1], [0]]) + assert TerminalExpr(grad(z), domain) == Matrix([[0], [0], [1]]) + assert TerminalExpr(grad(exp(1j*x)*sin(2*y)*z), domain) == Matrix([[1j*exp(1j*x)*sin(2*y)*z], [2*exp(1j*x)*cos(2*y)*z], [exp(1j*x)*sin(2*y)]]) + #============================================================================== @pytest.mark.skip(reason="New linearize() function does not accept 'LinearExpr' objects") def test_linearize_expr_2d_1(): diff --git a/sympde/exterior/tests/test_exterior.py b/sympde/exterior/tests/test_exterior.py index 07d28c88..d5778cb3 100644 --- a/sympde/exterior/tests/test_exterior.py +++ b/sympde/exterior/tests/test_exterior.py @@ -17,7 +17,7 @@ #============================================================================== def test_exterior_1(): - x, y, z = symbols('x y z') + x, y, z = symbols('x y z', real=True) a = Constant('a') n = Symbol('n') diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py index 07050c61..bce622f3 100644 --- a/sympde/topology/basic.py +++ b/sympde/topology/basic.py @@ -32,7 +32,7 @@ def coordinates(self): else: xyz = ['x', 'y', 'z'][:dim] - xyz = [Symbol(i) for i in xyz] + xyz = [Symbol(i, real=True) for i in xyz] self._coordinates = xyz if dim == 1: diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index f5d3537a..560024a6 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -104,7 +104,7 @@ def eval(cls, expr): return S.Zero elif isinstance(expr, Expr): - x = Symbol(cls.coordinate) + x = Symbol(cls.coordinate, real=True) if cls.logical: M = expr.atoms(Mapping) if len(M)>0: diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py index 1acf9a5d..31c852dd 100644 --- a/sympde/topology/domain.py +++ b/sympde/topology/domain.py @@ -839,7 +839,8 @@ def __new__(cls, name, dim, min_coords, max_coords): raise ValueError("Min coordinates must be smaller than max") coord_names = 'x1:{}'.format(dim + 1) - coordinates = symbols(coord_names) # [MCP 24.07.2024] we should be able to pass real=True but requires additional testing as it breaks some calls to TerminalExpr + + coordinates = symbols(coord_names, real=True) # Choose which type to use: # a) if dim <= 3, use Line, Square or Cube; diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 76090fc0..6efa5c90 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -184,7 +184,7 @@ def __new__(cls, name, dim=None, **kwargs): return obj if coordinates is None: - _coordinates = [Symbol(name) for name in ['x', 'y', 'z'][:pdim]] + _coordinates = [Symbol(name, real=True) for name in ['x', 'y', 'z'][:pdim]] else: if not isinstance(coordinates, (list, tuple, Tuple)): raise TypeError('> Expecting list, tuple, Tuple') @@ -193,7 +193,7 @@ def __new__(cls, name, dim=None, **kwargs): if not isinstance(a, (str, Symbol)): raise TypeError('> Expecting str or Symbol') - _coordinates = [Symbol(u) for u in coordinates] + _coordinates = [Symbol(u, real=True) for u in coordinates] obj._name = name obj._ldim = ldim @@ -203,16 +203,19 @@ def __new__(cls, name, dim=None, **kwargs): obj._is_minus = None obj._is_plus = None - lcoords = ['x1', 'x2', 'x3'][:ldim] - lcoords = [Symbol(i) for i in lcoords] - obj._logical_coordinates = Tuple(*lcoords) + lcoords_names = ['x1', 'x2', 'x3'][:ldim] + lcoords_symbols_real = [Symbol(i, real=True) for i in lcoords_names] + lcoords_symbols_real_dict = {key: symbol for key,symbol in zip(lcoords_names, lcoords_symbols_real)} + obj._logical_coordinates = Tuple(*lcoords_symbols_real) # coordinates must be symbols with real=True + lcoords_general_symbols = [Symbol(i) for i in lcoords_names] #list of symbols without real=True + # ... if not( obj._expressions is None ): - coords = ['x', 'y', 'z'][:pdim] + coords_names = ['x', 'y', 'z'][:pdim] # ... args = [] - for i in coords: + for i in coords_names: x = obj._expressions[i] x = sympify(x) args.append(x) @@ -225,13 +228,15 @@ def __new__(cls, name, dim=None, **kwargs): x = sympify(i) args = args.subs(x,0) # ... - - constants = list(set(args.free_symbols) - set(lcoords)) + # get constants by subtracting coordinates from list of free symbols + constants = list(set(args.free_symbols) - set(lcoords_general_symbols)) constants_values = {a.name:Constant(a.name) for a in constants} # subs constants as Constant objects instead of Symbol constants_values.update( kwargs ) d = {a:numpy_to_native_python(constants_values[a.name]) for a in constants} args = args.subs(d) + # subs coordinate symbols without real=True to symbols with real=True + args = args.subs(lcoords_symbols_real_dict) obj._expressions = args obj._constants = tuple(a for a in constants if isinstance(constants_values[a.name], Symbol)) @@ -919,7 +924,7 @@ def eval(cls, expr, domain, **options): if has(expr, DiffOperator): return cls( expr, domain, evaluate=False) else: - syms = symbols(ph_coords[:dim]) + syms = symbols(ph_coords[:dim], real=True) if isinstance(mapping, InterfaceMapping): mapping = mapping.minus # here we assume that the two mapped domains diff --git a/sympde/topology/tests/test_gallery.py b/sympde/topology/tests/test_gallery.py index 6222c9ce..5ce27784 100644 --- a/sympde/topology/tests/test_gallery.py +++ b/sympde/topology/tests/test_gallery.py @@ -4,7 +4,7 @@ from sympy import Tuple from sympy import symbols -x1, x2, x3 = symbols('x1, x2, x3') +x1, x2, x3 = symbols('x1, x2, x3', real=True) from sympde.topology import Interval, ProductDomain, InteriorDomain, Domain from sympde.topology import Line, Square, Cube, NCubeInterior diff --git a/sympde/topology/tests/test_logical_expr.py b/sympde/topology/tests/test_logical_expr.py index 99108c4d..59702e84 100644 --- a/sympde/topology/tests/test_logical_expr.py +++ b/sympde/topology/tests/test_logical_expr.py @@ -585,7 +585,7 @@ def test_symbolic_expr_3d_1(): def test_identity_mapping_2d_1(): dim = 2 - x1, x2 = symbols('x1, x2') + x1, x2 = symbols('x1, x2', real=True) M = IdentityMapping('M', dim=dim) @@ -626,7 +626,7 @@ def test_identity_mapping_2d_2(): def test_polar_mapping_2d_1(): dim = 2 - x1, x2 = symbols('x1, x2') + x1, x2 = symbols('x1, x2', real=True) constants = ['c1', 'c2', 'rmax', 'rmin'] c1, c2, rmax, rmin = [Constant(i) for i in constants] @@ -656,7 +656,7 @@ def test_polar_mapping_2d_1(): def test_target_mapping_2d_1(): dim = 2 - x1, x2 = symbols('x1, x2') + x1, x2 = symbols('x1, x2', real=True) constants = ['c1', 'c2', 'D', 'k'] c1, c2, D, k = [Constant(i) for i in constants] @@ -688,7 +688,7 @@ def test_target_mapping_2d_1(): def test_czarny_mapping_2d_1(): dim = 2 - x1, x2 = symbols('x1, x2') + x1, x2 = symbols('x1, x2', real=True) constants = ['c2', 'eps', 'b'] c2, eps, b = [Constant(i) for i in constants] @@ -725,7 +725,7 @@ def test_czarny_mapping_2d_1(): def test_collela_mapping_2d_1(): dim = 2 - x1, x2 = symbols('x1, x2') + x1, x2 = symbols('x1, x2', real=True) constants = ['eps', 'k1', 'k2'] eps, k1, k2 = [Constant(i) for i in constants] @@ -763,7 +763,7 @@ def test_collela_mapping_2d_1(): def test_torus_mapping_3d_1(): dim = 3 - x1, x2, x3 = symbols('x1, x2, x3') + x1, x2, x3 = symbols('x1, x2, x3', real=True) R0 = Constant('R0') M = TorusMapping('M', dim=dim) @@ -812,7 +812,7 @@ def test_torus_mapping_3d_1(): def test_twisted_target_mapping_3d_1(): dim = 3 - x1, x2, x3 = symbols('x1, x2, x3') + x1, x2, x3 = symbols('x1, x2, x3', real=True) constants = ['c1', 'c2', 'c3', 'D', 'k'] c1, c2, c3, D, k = [Constant(i) for i in constants]