diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 61fb8c34..f06b7ef3 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -84,3 +84,37 @@ jobs: working-directory: ./scratch run: | coverage report --ignore-errors --show-missing --sort=cover + + minimum-sympy: + name: Ubuntu / Python 3.11 / SymPy 1.12 + runs-on: ubuntu-24.04 + + steps: + - uses: actions/checkout@v6 + + - name: Set up Python 3.11 + uses: actions/setup-python@v6 + with: + python-version: '3.11' + cache: 'pip' + cache-dependency-path: | + pyproject.toml + + - name: Upgrade pip, setuptools, and wheel + run: | + pip install --upgrade pip setuptools wheel + + - name: Install project with minimum SymPy + run: | + pip install .[test] 'sympy == 1.12' + pip freeze + + - name: Initialize test directory + run: | + mkdir scratch + + - name: Run tests + working-directory: ./scratch + run: >- + pytest -n auto --dist loadgroup + --pyargs sympde -ra diff --git a/CHANGELOG.md b/CHANGELOG.md index 89430b40..dd8c8c7d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,13 +14,19 @@ All notable changes to this project will be documented in this file. ### Changed - Convert the README to Markdown and expand the installation, testing, and documentation-build instructions. +- #106: Require SymPy 1.12 or newer and support maintained SymPy releases below 1.15. +- #106: Update SymPDE's symbolic objects for modern SymPy's structural arguments, traversal, substitution, and reconstruction requirements. - [DEVELOPER] Modernize the package metadata and consolidate the pytest and coverage configuration in `pyproject.toml`. - [DEVELOPER] Modernize the Sphinx configuration, API-documentation generation, and API docstring markup. +- [DEVELOPER] Test the minimum supported SymPy version in a dedicated CI job. ### Fixed - #140: Fix the rendering of the `Domain.join()` API documentation. - #166: Restore warning-free documentation builds, including bibliography and mathematical notation support. +- #106: Preserve symbol assumptions and free-symbol information with modern SymPy. +- #106: Keep same-named functions from different function spaces distinct during symbolic simplification. +- #106: Fix form-linearity checks involving integral sums with modern SymPy. ### Removed diff --git a/pyproject.toml b/pyproject.toml index d58959e4..bbb33757 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ classifiers = [ "Topic :: Scientific/Engineering :: Mathematics", ] dependencies = [ - 'sympy >= 1.5, < 1.10', + 'sympy >= 1.12, < 1.15', 'h5py', 'pyyaml', 'numpy', diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index ac59ef41..7d40dab3 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -102,8 +102,8 @@ from sympy.core.containers import Tuple from sympy.core.singleton import S from sympy.core.decorators import call_highest_priority +from sympy.utilities.iterables import is_sequence -from sympde.old_sympy_utilities import is_sequence from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery from sympde.topology.space import ScalarFunction, VectorFunction @@ -1860,4 +1860,3 @@ def add_basicop(expr): Basic._constructor_postprocessor_mapping[BasicOperator] = { "Add": [add_basicop], } - diff --git a/sympde/calculus/matrices.py b/sympde/calculus/matrices.py index 16b827e1..c796e638 100644 --- a/sympde/calculus/matrices.py +++ b/sympde/calculus/matrices.py @@ -2,7 +2,7 @@ from sympy import Add, Mul, Pow from sympy import sympify from sympy.core.decorators import call_highest_priority -from sympde.core.basic import _coeffs_registery, Basic +from sympde.core.basic import _coeffs_registery, _new_basic, Basic class MatrixSymbolicExpr(Expr): is_commutative = False @@ -275,7 +275,7 @@ def _sympystr(self, printer): class MatrixElement(Expr): def __new__(cls, base, indices, **options): - return Expr.__new__(cls, base, indices) + return _new_basic(cls, base, indices) @property def base(self): @@ -293,4 +293,3 @@ def _sympystr(self, printer): "Mul": [lambda x: MatSymbolicMul(*x.args)], "Add": [lambda x: MatSymbolicAdd(*x.args)] } - diff --git a/sympde/core/algebra.py b/sympde/core/algebra.py index 7b626117..acebf082 100644 --- a/sympde/core/algebra.py +++ b/sympde/core/algebra.py @@ -33,8 +33,8 @@ from sympy import S from sympy import Basic from sympy import Indexed, IndexedBase +from sympy.utilities.iterables import is_sequence -from sympde.old_sympy_utilities import is_sequence from .basic import CalculusFunction from .basic import _coeffs_registery diff --git a/sympde/core/basic.py b/sympde/core/basic.py index 5f8fff14..c7569ac8 100644 --- a/sympde/core/basic.py +++ b/sympde/core/basic.py @@ -3,9 +3,78 @@ from sympy import Function from sympy import Number from sympy import NumberSymbol +from sympy import Integer +from sympy import sympify from sympy.core import Basic from sympy.core import Symbol -from sympy.tensor import IndexedBase +from sympy.core.containers import Dict, Tuple +from sympy.core.symbol import Str + + +class _NoneArgument(Basic): + """SymPy-compatible representation of ``None`` in expression arguments.""" + + +_none = _NoneArgument() + + +def _sympify_argument(value): + """Convert structural metadata to objects accepted in ``Basic.args``.""" + if isinstance(value, Basic): + return value + if value is None: + return _none + if isinstance(value, str): + return Str(value) + if isinstance(value, dict): + return Dict(*((_sympify_argument(k), _sympify_argument(v)) + for k, v in value.items())) + if isinstance(value, (tuple, list)): + return Tuple(*(_sympify_argument(v) for v in value)) + return sympify(value) + + +def _new_basic(cls, *args, **options): + """Construct a ``Basic`` after normalizing all structural arguments.""" + args = tuple(_sympify_argument(arg) for arg in args) + return Basic.__new__(cls, *args, **options) + + +def _is_none_argument(value): + return value is _none + + +def _restore_optional_argument(value): + """Restore an optional value stored in a symbolic argument list.""" + return None if _is_none_argument(value) else value + + +def _restore_string_argument(value): + """Restore a Python string stored as a SymPy ``Str`` object.""" + return value.name if isinstance(value, Str) else value + + +def _restore_int_argument(value): + """Restore a Python integer stored as a SymPy ``Integer`` object.""" + value = _restore_optional_argument(value) + return int(value) if isinstance(value, Integer) else value + + +def _restore_container_argument(value): + """Restore Python containers used as non-expression metadata.""" + value = _restore_optional_argument(value) + if isinstance(value, Str): + return value.name + if isinstance(value, Dict): + return { + _restore_container_argument(k): _restore_container_argument(v) + for k, v in value.items() + } + if isinstance(value, Tuple): + return tuple(_restore_container_argument(v) for v in value) + if isinstance(value, Integer): + return int(value) + return value #============================================================================== class Constant(Symbol): @@ -34,13 +103,6 @@ class CalculusFunction(Function): functions when manipulating our expressions""" pass -#============================================================================== -class BasicMapping(IndexedBase): - """ - Represents a basic class for mapping. - """ - pass - #============================================================================== class BasicDerivable(Basic): pass diff --git a/sympde/expr/equation.py b/sympde/expr/equation.py index 5b4ed1aa..0a92b71d 100644 --- a/sympde/expr/equation.py +++ b/sympde/expr/equation.py @@ -13,6 +13,11 @@ from sympde.topology import Trace, trace_0, trace_1 from sympde.calculus import grad, dot from sympde.core.utils import random_string +from sympde.core.basic import ( + _is_none_argument, + _new_basic, + _restore_optional_argument, +) from .expr import BilinearForm, LinearForm from .expr import linearize @@ -191,6 +196,8 @@ def rhs(self): class Equation(Basic): def __new__(cls, lhs, rhs, trials, tests, bc=None, constraint=None): + bc = _restore_optional_argument(bc) + constraint = _restore_optional_argument(constraint) # ... if not isinstance(lhs, BilinearForm): raise UnconsistentLhsError('> lhs must be a bilinear') @@ -332,7 +339,7 @@ def __new__(cls, lhs, rhs, trials, tests, bc=None, constraint=None): tests = Tuple(*tests) # ... - return Basic.__new__(cls, lhs, rhs, trials, tests, bc, constraint) + return _new_basic(cls, lhs, rhs, trials, tests, bc, constraint) @property def lhs(self): @@ -352,11 +359,13 @@ def test_functions(self): @property def bc(self): - return self._args[4] + bc = self._args[4] + return None if _is_none_argument(bc) else bc @property def constraint(self): - return self._args[5] + constraint = self._args[5] + return None if _is_none_argument(constraint) else constraint #============================================================================== # TODO must subtitute expr by given args => call then create BasicForm diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index aa148b04..d6e026d4 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -8,6 +8,7 @@ from sympde.core.basic import _coeffs_registery from sympde.core.basic import CalculusFunction +from sympde.core.basic import _new_basic from sympde.core.algebra import (Dot_1d, Dot_2d, Inner_2d, Cross_2d, Dot_3d, Inner_3d, Cross_3d) @@ -973,15 +974,15 @@ def _tensorize_atomic_expr(expr, d_atoms): class Basic1dForm(AtomicExpr): def __new__(cls, name, axis, weight=S.One): - return Basic.__new__(cls, name, axis, weight) + return _new_basic(cls, name, axis, weight) @property def name(self): - return self._args[0] + return self._args[0].name @property def axis(self): - return self._args[1] + return int(self._args[1]) @property def weight(self): @@ -1003,38 +1004,48 @@ def _sympystr(self, printer): class Mass(Basic1dForm): - def __new__(cls, axis, weight=S.One): + def __new__(cls, axis, weight=S.One, *structural_args): + if structural_args: + return Basic1dForm.__new__(cls, axis, weight, *structural_args) # name = 'Mass' name = 'M' - return Basic.__new__(cls, name, axis, weight) + return Basic1dForm.__new__(cls, name, axis, weight) class Stiffness(Basic1dForm): - def __new__(cls, axis, weight=S.One): + def __new__(cls, axis, weight=S.One, *structural_args): + if structural_args: + return Basic1dForm.__new__(cls, axis, weight, *structural_args) # name = 'Stiffness' name = 'S' - return Basic.__new__(cls, name, axis, weight) + return Basic1dForm.__new__(cls, name, axis, weight) class Advection(Basic1dForm): - def __new__(cls, axis, weight=S.One): + def __new__(cls, axis, weight=S.One, *structural_args): + if structural_args: + return Basic1dForm.__new__(cls, axis, weight, *structural_args) # name = 'Advection' name = 'A' - return Basic.__new__(cls, name, axis, weight) + return Basic1dForm.__new__(cls, name, axis, weight) class AdvectionT(Basic1dForm): - def __new__(cls, axis, weight=S.One): + def __new__(cls, axis, weight=S.One, *structural_args): + if structural_args: + return Basic1dForm.__new__(cls, axis, weight, *structural_args) # name = 'AdvectionT' name = 'AT' - return Basic.__new__(cls, name, axis, weight) + return Basic1dForm.__new__(cls, name, axis, weight) class Bilaplacian(Basic1dForm): - def __new__(cls, axis, weight=S.One): + def __new__(cls, axis, weight=S.One, *structural_args): + if structural_args: + return Basic1dForm.__new__(cls, axis, weight, *structural_args) # name = 'Bilaplacian' name = 'B' - return Basic.__new__(cls, name, axis, weight) + return Basic1dForm.__new__(cls, name, axis, weight) Mass_0 = Mass(0) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 67557467..a06acd58 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -9,8 +9,8 @@ from sympy.core import Expr, Add, Mul from sympy.core.numbers import Zero as sy_Zero from sympy.core.containers import Tuple +from sympy.utilities.iterables import is_sequence -from sympde.old_sympy_utilities import is_sequence from sympde.core.basic import CalculusFunction from sympde.core.basic import Constant from sympde.core.utils import random_string @@ -402,7 +402,7 @@ def __call__(self, *tests, **kwargs): # Make sure that 'values' is always a list if len(tests) == 1: values = tests[0] - if not is_sequence(values, vector=isinstance(values, VectorFunction)): + if isinstance(values, VectorFunction) or not is_sequence(values): values = [values] else: values = tests @@ -513,8 +513,8 @@ def __call__(self, trials, tests, **kwargs): expr = self._update_free_variables(**kwargs) # If needed, convert positional arguments to lists - if not is_sequence(trials, vector=isinstance(trials, VectorFunction)): trials = [trials] - if not is_sequence(tests, vector=isinstance(tests, VectorFunction)): tests = [tests ] + if isinstance(trials, VectorFunction) or not is_sequence(trials): trials = [trials] + if isinstance(tests, VectorFunction) or not is_sequence(tests): tests = [tests ] # Concatenate input values into single list values = [*trials, *tests] @@ -855,7 +855,7 @@ def mul_integral(expr): def mul_add(expr): args_int = [a for a in expr.args if isinstance(a, IntAdd)] coeff = Mul(*[a for a in expr.args if not a in args_int]) - args_int[0] = [IntAdd(*[coeff*i for i in args_int[0].args])] + args_int[0] = IntAdd(*[coeff*i for i in args_int[0].args]) return Mul(*args_int) def add_int(expr): diff --git a/sympde/expr/tests/test_equation_2d.py b/sympde/expr/tests/test_equation_2d.py index cb446f2f..fa82ff26 100644 --- a/sympde/expr/tests/test_equation_2d.py +++ b/sympde/expr/tests/test_equation_2d.py @@ -109,6 +109,7 @@ def test_equation_2d_1(): a = BilinearForm((v,u), a1(v,u) + a_B1(v,u)) l = LinearForm(v, l1(v) + l2(v)) equation = Equation(a, l, tests=v, trials=u) + assert equation.func(*equation.args) == equation # ... # ... @@ -415,4 +416,3 @@ def teardown_module(): def teardown_function(): from sympy.core import cache cache.clear_cache() - diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 77a96008..966c2223 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -14,7 +14,7 @@ from sympde.topology import dx1, dx2, dx3 from sympde.topology import dx, dy, dz -from sympde.topology import Mapping +from sympde.topology import SymbolicMapping from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace from sympde.topology import element_of, elements_of from sympde.topology import InteriorDomain, Union @@ -658,7 +658,7 @@ def test_terminal_expr_bilinear_2d_4(): def test_terminal_expr_bilinear_3d_1(): domain = Domain('Omega', dim=3) - M = Mapping('M', dim=3) + M = SymbolicMapping('M', dim=3) mapped_domain = M(domain) @@ -687,6 +687,72 @@ 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_kernel_expression_structure_linear_2d(): + """Check that TerminalExpr splits a mixed interior/boundary LinearForm into + the expected number and kind of KernelExpression objects.""" + + domain = Domain('Omega', dim=2) + B1 = Boundary(r'\Gamma_1', domain) + + x, y = domain.coordinates + nn = NormalVector('nn') + + V = ScalarFunctionSpace('V', domain) + v = element_of(V, name='v') + + int_0 = lambda expr: integral(domain, expr) + int_1 = lambda expr: integral(B1, expr) + + g = Matrix((x**2, y**2)) + l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v)) + + kernel_expr = TerminalExpr(l, domain) + assert len(kernel_expr) == 2 + + domain_exprs = [k for k in kernel_expr if isinstance(k.target, InteriorDomain)] + boundary_exprs = [k for k in kernel_expr if isinstance(k.target, Boundary)] + assert len(domain_exprs) == 1 + assert len(boundary_exprs) == 1 + + assert domain_exprs[0].target == domain.interior + assert boundary_exprs[0].target == B1 + + assert domain_exprs[0].expr == x*y*v + assert boundary_exprs[0].expr == v*(x**2*nn[0] + y**2*nn[1]) + +#============================================================================== +def test_terminal_expr_kernel_expression_structure_bilinear_2d(): + """Check that TerminalExpr splits a mixed interior/boundary BilinearForm + into the expected number and kind of KernelExpression objects.""" + + domain = Domain('Omega', dim=2) + B1 = Boundary(r'\Gamma_1', domain) + + nn = NormalVector('nn') + + V = ScalarFunctionSpace('V', domain) + u, v = elements_of(V, names='u, v') + + int_0 = lambda expr: integral(domain, expr) + int_1 = lambda expr: integral(B1, expr) + + a = BilinearForm((u, v), int_0(u*v + dot(grad(u), grad(v))) + int_1(v*dot(grad(u), nn))) + + kernel_expr = TerminalExpr(a, domain) + assert len(kernel_expr) == 2 + + domain_exprs = [k for k in kernel_expr if isinstance(k.target, InteriorDomain)] + boundary_exprs = [k for k in kernel_expr if isinstance(k.target, Boundary)] + assert len(domain_exprs) == 1 + assert len(boundary_exprs) == 1 + + assert domain_exprs[0].target == domain.interior + assert boundary_exprs[0].target == B1 + + assert domain_exprs[0].expr == u*v + dx1(u)*dx1(v) + dx2(u)*dx2(v) + assert boundary_exprs[0].expr == v*(dx1(u)*nn[0] + dx2(u)*nn[1]) + #============================================================================== def test_terminal_expr_1d_1(): domain = Domain('Omega', dim=1) diff --git a/sympde/expr/tests/test_tensor_2d.py b/sympde/expr/tests/test_tensor_2d.py index 0fc9e31d..7a68e429 100644 --- a/sympde/expr/tests/test_tensor_2d.py +++ b/sympde/expr/tests/test_tensor_2d.py @@ -12,7 +12,7 @@ from sympde.topology import elements_of from sympde.topology import Boundary from sympde.topology import Domain -from sympde.topology import Mapping +from sympde.topology import SymbolicMapping from sympde.expr.expr import BilinearForm, integral from sympde.expr.evaluation import TensorExpr @@ -63,7 +63,7 @@ def test_tensorize_2d_1_mapping(): DIM = 2 - M = Mapping('Map', dim=DIM) + M = SymbolicMapping('Map', dim=DIM) domain = M(Domain('Omega', dim=DIM)) @@ -88,7 +88,7 @@ def test_tensorize_2d_1_mapping(): def test_tensorize_2d_2_mapping(): DIM = 2 - M = Mapping('M', dim=DIM) + M = SymbolicMapping('M', dim=DIM) domain = M(Domain('Omega', dim=DIM)) V = VectorFunctionSpace('V', domain) diff --git a/sympde/exterior/calculus.py b/sympde/exterior/calculus.py index c9b97059..d3bf6bf9 100644 --- a/sympde/exterior/calculus.py +++ b/sympde/exterior/calculus.py @@ -14,6 +14,7 @@ from sympy.core.singleton import Singleton from sympy.core import Add, Mul from sympy.core.singleton import S +from sympy.utilities.iterables import is_sequence from sympde.core.basic import _coeffs_registery from sympde.core import LinearOperator diff --git a/sympde/exterior/compiler.py b/sympde/exterior/compiler.py index 0fefba08..c67c6ad5 100644 --- a/sympde/exterior/compiler.py +++ b/sympde/exterior/compiler.py @@ -7,6 +7,7 @@ from sympy import S from sympy import Indexed from sympy import Tuple +from sympy.utilities.iterables import is_sequence from sympde.core.basic import _coeffs_registery from sympde.core.basic import CalculusFunction diff --git a/sympde/exterior/datatype.py b/sympde/exterior/datatype.py index cb96ba5e..09f2f033 100644 --- a/sympde/exterior/datatype.py +++ b/sympde/exterior/datatype.py @@ -14,12 +14,11 @@ from sympy.core import Add, Mul from sympy.core.singleton import S -from sympde.old_sympy_utilities import with_metaclass from sympde.core.basic import _coeffs_registery from sympde.core import LinearOperator #============================================================================== -class FormType(with_metaclass(Singleton, Basic)): +class FormType(Basic, metaclass=Singleton): """Base class representing differential form types""" _index = None diff --git a/sympde/exterior/form.py b/sympde/exterior/form.py index dc4652e7..b45b9f76 100644 --- a/sympde/exterior/form.py +++ b/sympde/exterior/form.py @@ -34,16 +34,18 @@ def __new__(cls, name, index, dim): index = get_index_form(index) - return Basic.__new__(cls, name, index, dim) - - @property - def name(self): - return self._args[0] + obj = Symbol.__xnew__(cls, name) + obj._index = index + obj._dim = Integer(dim) if isinstance(dim, int) else dim + return obj @property def index(self): - return self._args[1] + return self._index @property def dim(self): - return self._args[2] + return int(self._dim) if self._dim.is_Integer else self._dim + + def _hashable_content(self): + return Symbol._hashable_content(self) + (self.index, self._dim) diff --git a/sympde/exterior/tests/test_exterior.py b/sympde/exterior/tests/test_exterior.py index d5778cb3..944c9e8b 100644 --- a/sympde/exterior/tests/test_exterior.py +++ b/sympde/exterior/tests/test_exterior.py @@ -13,6 +13,17 @@ from sympde.exterior import DifferentialForm +def test_differential_form_symbol_identity(): + u = DifferentialForm('u', index=1, dim=3) + + assert u == DifferentialForm('u', index=1, dim=3) + assert hash(u) == hash(DifferentialForm('u', index=1, dim=3)) + assert u != DifferentialForm('u', index=2, dim=3) + assert u != DifferentialForm('u', index=1, dim=2) + assert u.free_symbols == {u} + assert u.args == () + + #============================================================================== def test_exterior_1(): diff --git a/sympde/printing/latex.py b/sympde/printing/latex.py index 28c89770..50c7b199 100644 --- a/sympde/printing/latex.py +++ b/sympde/printing/latex.py @@ -202,6 +202,9 @@ def _print_Kron(self, expr): def _print_Nil(self, expr): return r'\ldots' + def _print_BasicDomain(self, expr): + return translate(expr.name) + def _print_Domain(self, expr): return translate(expr.name) diff --git a/sympde/printing/tests/test_latex.py b/sympde/printing/tests/test_latex.py index 1c8cef5c..c06a9fe8 100644 --- a/sympde/printing/tests/test_latex.py +++ b/sympde/printing/tests/test_latex.py @@ -4,7 +4,7 @@ from sympde.calculus import grad, dot, inner, rot, div #from sympde.topology import (dx, dy, dz) -from sympde.topology import Domain, Boundary +from sympde.topology import Domain, Boundary, Square from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace from sympde.topology import element_of from sympde.expr import BilinearForm, LinearForm, integral @@ -13,6 +13,19 @@ from sympde.printing.latex import latex +def test_latex_domain_subclasses_use_their_names(): + patch_a = Square('A') + patch_b = Square('B') + domain = Domain.join( + [patch_a, patch_b], + [((0, 0, 1), (1, 0, -1), 1)], + 'Omega', + ) + + assert latex(patch_a.interior) == 'A' + assert latex(domain.interfaces) == 'A|B' + + #============================================================================== def test_latex_1d(): @@ -300,4 +313,3 @@ def teardown_module(): def teardown_function(): from sympy.core import cache cache.clear_cache() - diff --git a/sympde/topology/analytical_mapping.py b/sympde/topology/analytical_mapping.py index 05b9760b..c5097a30 100644 --- a/sympde/topology/analytical_mapping.py +++ b/sympde/topology/analytical_mapping.py @@ -1,6 +1,6 @@ -from .mapping import Mapping +from .mapping import AnalyticMapping -class IdentityMapping(Mapping): +class IdentityMapping(AnalyticMapping): """ Represents an identity 1D/2D/3D Mapping object. @@ -12,7 +12,7 @@ class IdentityMapping(Mapping): 'z': 'x3'} #============================================================================== -class AffineMapping(Mapping): +class AffineMapping(AnalyticMapping): """ Represents a 1D/2D/3D Affine Mapping object. @@ -24,7 +24,7 @@ class AffineMapping(Mapping): 'z': 'c3 + a31*x1 + a32*x2 + a33*x3'} #============================================================================== -class PolarMapping(Mapping): +class PolarMapping(AnalyticMapping): """ Represents a Polar 2D Mapping object (Annulus). @@ -38,7 +38,7 @@ class PolarMapping(Mapping): _pdim = 2 #============================================================================== -class TargetMapping(Mapping): +class TargetMapping(AnalyticMapping): """ Represents a Target 2D Mapping object. @@ -52,7 +52,7 @@ class TargetMapping(Mapping): _pdim = 2 #============================================================================== -class CzarnyMapping(Mapping): +class CzarnyMapping(AnalyticMapping): """ Represents a Czarny 2D Mapping object. @@ -67,7 +67,7 @@ class CzarnyMapping(Mapping): _pdim = 2 #============================================================================== -class CollelaMapping2D(Mapping): +class CollelaMapping2D(AnalyticMapping): """ Represents a Collela 2D Mapping object. @@ -79,7 +79,7 @@ class CollelaMapping2D(Mapping): _pdim = 2 #============================================================================== -class TorusMapping(Mapping): +class TorusMapping(AnalyticMapping): """ Parametrization of a torus (or a portion of it) of major radius R0, using toroidal coordinates (x1, x2, x3) = (r, theta, phi), where: @@ -98,7 +98,7 @@ class TorusMapping(Mapping): #============================================================================== # TODO [YG, 07.10.2022]: add test in sympde/topology/tests/test_logical_expr.py -class TorusSurfaceMapping(Mapping): +class TorusSurfaceMapping(AnalyticMapping): """ 3D surface obtained by "slicing" the torus above at r = a. The parametrization uses the coordinates (x1, x2) = (theta, phi), where: @@ -116,7 +116,7 @@ class TorusSurfaceMapping(Mapping): #============================================================================== # TODO [YG, 07.10.2022]: add test in sympde/topology/tests/test_logical_expr.py -class TwistedTargetSurfaceMapping(Mapping): +class TwistedTargetSurfaceMapping(AnalyticMapping): """ 3D surface obtained by "twisting" the TargetMapping out of the (x, y) plane @@ -129,7 +129,7 @@ class TwistedTargetSurfaceMapping(Mapping): _pdim = 3 #============================================================================== -class TwistedTargetMapping(Mapping): +class TwistedTargetMapping(AnalyticMapping): """ 3D volume obtained by "extruding" the TwistedTargetSurfaceMapping along z. @@ -142,7 +142,7 @@ class TwistedTargetMapping(Mapping): _pdim = 3 #============================================================================== -class SphericalMapping(Mapping): +class SphericalMapping(AnalyticMapping): """ Parametrization of a sphere (or a portion of it) using spherical coordinates (x1, x2, x3) = (r, theta, phi), where: diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py index 4eb371d6..8bee0976 100644 --- a/sympde/topology/basic.py +++ b/sympde/topology/basic.py @@ -8,6 +8,15 @@ from sympy.core.containers import Tuple from sympy.tensor import IndexedBase +from sympde.core.basic import ( + _is_none_argument, + _new_basic, + _restore_container_argument, + _restore_int_argument, + _restore_optional_argument, + _restore_string_argument, +) + #============================================================================== class BasicDomain(Basic): _dim = None @@ -52,8 +61,15 @@ class InteriorDomain(BasicDomain): Examples """ - def __new__(cls, name, dim=None, dtype=None, mapping=None, logical_domain=None): - target = None + def __new__(cls, name, dim=None, dtype=None, mapping=None, + logical_domain=None, target=None): + name = _restore_string_argument(name) + dim = _restore_int_argument(dim) + dtype = _restore_container_argument(dtype) + mapping = _restore_optional_argument(mapping) + logical_domain = _restore_optional_argument(logical_domain) + target = _restore_optional_argument(target) + if not isinstance(name, str): target = name name = name.name @@ -64,7 +80,8 @@ def __new__(cls, name, dim=None, dtype=None, mapping=None, logical_domain=None): assert mapping is None and logical_domain is None or \ mapping is not None and logical_domain is not None - obj = Basic.__new__(cls, name) + obj = _new_basic(cls, name, dim, dtype, mapping, + logical_domain, target) obj._dim = dim obj._target = target @@ -76,7 +93,7 @@ def __new__(cls, name, dim=None, dtype=None, mapping=None, logical_domain=None): @property def name(self): - return self.args[0] + return self.args[0].name @property def target(self): @@ -231,13 +248,17 @@ class Interval(InteriorDomain): _dim = 1 def __new__(cls, name=None, coordinate=None, bounds=None): + name = _restore_string_argument(name) + coordinate = _restore_optional_argument(coordinate) + bounds = _restore_container_argument(bounds) + if name is None: name = 'Interval' if bounds is None: bounds = (0, 1) - obj = Basic.__new__(cls, name) + obj = _new_basic(cls, name, coordinate, bounds) if coordinate: obj._coordinates = [coordinate] @@ -247,7 +268,7 @@ def __new__(cls, name=None, coordinate=None, bounds=None): @property def name(self): - return self.args[0] + return self.args[0].name @property def bounds(self): @@ -264,13 +285,20 @@ class Boundary(BasicDomain): """ def __new__(cls, name, domain, axis=None, ext=None, mapping=None, logical_domain=None): + name = _restore_string_argument(name) + axis = _restore_int_argument(axis) + ext = _restore_int_argument(ext) + mapping = _restore_optional_argument(mapping) + logical_domain = _restore_optional_argument(logical_domain) + if axis is not None: assert isinstance(axis, int) if ext is not None: assert isinstance(ext, int) - obj = Basic.__new__(cls, name, domain, axis, ext) + obj = _new_basic(cls, name, domain, axis, ext, + mapping, logical_domain) obj._mapping = mapping obj._logical_domain = logical_domain @@ -278,7 +306,7 @@ def __new__(cls, name, domain, axis=None, ext=None, mapping=None, logical_domain @property def name(self): - return self.args[0] + return self.args[0].name @property def domain(self): @@ -286,11 +314,13 @@ def domain(self): @property def axis(self): - return self.args[2] + axis = self.args[2] + return None if _is_none_argument(axis) else int(axis) @property def ext(self): - return self.args[3] + ext = self.args[3] + return None if _is_none_argument(ext) else int(ext) @property def mapping(self): @@ -450,7 +480,7 @@ class Interface(BasicDomain): bnd_plus : Boundary Boundary on the "plus" side of the interface. - mapping : Mapping, optional + mapping : SymbolicMapping, optional Mapping from the logical domain to the physical domain, if available. logical_domain : BasicDomain, optional @@ -470,7 +500,13 @@ class Interface(BasicDomain): and T. Dokken, E. Quak, V. Skytt. Requirements from Isogeometric Analysis for changes in product design ontologies, 2010. """ - def __new__(cls, name, bnd_minus, bnd_plus, *, mapping=None, logical_domain=None, ornt=None): + def __new__(cls, name, bnd_minus, bnd_plus, ornt=None, + mapping=None, logical_domain=None): + + name = _restore_string_argument(name) + ornt = _restore_container_argument(ornt) + mapping = _restore_optional_argument(mapping) + logical_domain = _restore_optional_argument(logical_domain) if not isinstance(name , str ): raise TypeError(name) if not isinstance(bnd_minus, Boundary): raise TypeError(bnd_minus) @@ -489,9 +525,12 @@ def __new__(cls, name, bnd_minus, bnd_plus, *, mapping=None, logical_domain=None # If provided, check that mapping is consistent with the boundaries if mapping is not None: - from sympde.topology.mapping import Mapping - if not isinstance(mapping, Mapping): - raise TypeError(f'mapping must be of type Mapping, got {type(mapping)} instead') + # SymbolicMapping, not Mapping: an interface can carry an + # InterfaceMapping (a StructuralMapping) here, which since WP06d-4a + # no longer subclasses Mapping. + from sympde.topology.mapping import SymbolicMapping + if not isinstance(mapping, SymbolicMapping): + raise TypeError(f'mapping must be a SymbolicMapping, got {type(mapping)} instead') if mapping.ldim != ldim: raise ValueError(f'mapping should have logical dimension = {ldim}, got {mapping.ldim} instead') @@ -518,7 +557,8 @@ def __new__(cls, name, bnd_minus, bnd_plus, *, mapping=None, logical_domain=None # TODO [YG 10.02.2026]: relax this requirement ASAP assert bnd_minus.axis == bnd_plus.axis - obj = Basic.__new__(cls, name, bnd_minus, bnd_plus, ornt) + obj = _new_basic(cls, name, bnd_minus, bnd_plus, ornt, + mapping, logical_domain) obj._mapping = mapping obj._logical_domain = logical_domain return obj @@ -529,7 +569,7 @@ def dim(self): @property def name(self): - return self.args[0] + return self.args[0].name @property def minus(self): @@ -541,7 +581,12 @@ def plus(self): @property def ornt(self): - return self.args[3] + ornt = self.args[3] + if _is_none_argument(ornt): + return None + if isinstance(ornt, Tuple): + return tuple(int(o) for o in ornt) + return int(ornt) @property def axis(self): @@ -620,6 +665,7 @@ def __setitem__(self, key, value): self._data[key] = value # ========================================== + # abstract methods # ========================================== def __getitem__(self, key): @@ -640,4 +686,3 @@ def __lt__(self, other): # ========================================== - diff --git a/sympde/topology/callable_mapping.py b/sympde/topology/callable_mapping.py index f1f4d60f..1768214d 100644 --- a/sympde/topology/callable_mapping.py +++ b/sympde/topology/callable_mapping.py @@ -1,95 +1,20 @@ -from sympy import Symbol - -from sympde.utilities.utils import lambdify_sympde -from .mapping import Mapping, BasicCallableMapping - -__all__ = ('CallableMapping',) - -#============================================================================== -class CallableMapping(BasicCallableMapping): - - def __init__( self, mapping, **kwargs ): - - # Extract information from class - assert isinstance(mapping, Mapping) - - variables = mapping.logical_coordinates - expressions = mapping.expressions - - jac = mapping.jacobian_expr - inv_jac = mapping.jacobian_inv_expr - - metric = mapping.metric_expr - metric_det = mapping.metric_det_expr - - constants = mapping.constants - params = {a.name:a for a in constants} - params.update( kwargs ) - for p in params.values(): - assert not isinstance(p, Symbol) - - if params: - subs = {a:params[a.name] for a in constants} - expressions = expressions.subs(subs) - jac = jac.subs(subs) - inv_jac = inv_jac.subs(subs) - metric = metric.subs(subs) - metric_det = metric_det.subs(subs) - - # Callable function: __call__ - self._func_eval = tuple(lambdify_sympde( variables, expr) for expr in expressions) - - # Callable function: jac_mat - self._jacobian = lambdify_sympde( variables, jac) - - # Callable function: jac_mat_inv - self._jacobian_inv = lambdify_sympde( variables, inv_jac) - - # Callable function: metric - self._metric = lambdify_sympde( variables, metric) - - # Callable function: metric_det - self._metric_det = lambdify_sympde( variables, metric_det) - - # Symbolic information - self._params = params - self._symbolic_mapping = mapping - - #-------------------------------------------------------------------------- - # Abstract interface - #-------------------------------------------------------------------------- - def __call__(self, *eta): - return tuple( f( *eta ) for f in self._func_eval) - - def jacobian(self, *eta): - return self._jacobian( *eta ) - - def jacobian_inv(self, *eta): - """ Compute the inverse Jacobian matrix, if possible.""" - return self._jacobian_inv(*eta) - - def metric(self, *eta): - return self._metric( *eta ) - - def metric_det(self, *eta): - return self._metric_det( *eta ) - - @property - def ldim(self): - return self.symbolic_mapping.ldim - - @property - def pdim(self): - return self.symbolic_mapping.pdim - - #-------------------------------------------------------------------------- - # Symbolic information - #-------------------------------------------------------------------------- - @property - def params( self ): - return self._params - - @property - def symbolic_mapping( self ): - return self._symbolic_mapping - +# The point-evaluable analytic wrapper ``CallableMapping`` was removed in +# WP06d-2: an ``AnalyticMapping`` lambdifies its own expressions (WP06c) and is +# its own callable mapping. This module now only re-exports +# ``BasicCallableMapping`` (the plain-Python point-evaluation ABC) for the +# psydac call sites that still import it from here; those migrate to +# ``DefinedMapping`` in WP06d-3. + +from .mapping import BasicCallableMapping + +__all__ = ('BasicCallableMapping',) + + +def __getattr__(name): + if name == 'CallableMapping': + raise AttributeError( + "CallableMapping was removed in sympde WP06d-2. Analytic mappings " + "are point-evaluable directly: subclass AnalyticMapping (an " + "AnalyticMapping instance is its own callable mapping), or wrap a " + "callable with DiscreteMapping(F, name).") + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/sympde/topology/datatype.py b/sympde/topology/datatype.py index 6ba5d4b2..3d7722bf 100644 --- a/sympde/topology/datatype.py +++ b/sympde/topology/datatype.py @@ -13,10 +13,9 @@ from sympy.core import Expr from sympy.core.containers import Tuple from sympy.core.singleton import Singleton -from sympde.old_sympy_utilities import with_metaclass #============================================================================== -class SpaceType(with_metaclass(Singleton, Basic)): +class SpaceType(Basic, metaclass=Singleton): """Base class representing function space types""" pass @@ -48,7 +47,7 @@ class UndefinedSpaceType(SpaceType): 'undefined': UndefinedSpace} #============================================================================== -class RegularityType(with_metaclass(Singleton, Basic)): +class RegularityType(Basic, metaclass=Singleton): """Base class representing the regularity of a space of functions""" _index = None diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index 560024a6..84a0a73c 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -18,11 +18,10 @@ from sympy import log from sympy import preorder_traversal from sympy import cacheit -from sympde.old_sympy_utilities import is_sequence +from sympy.utilities.iterables import is_sequence from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery -from sympde.core.basic import BasicMapping from sympde.core.algebra import LinearOperator from sympde.calculus.core import minus, plus from sympde.calculus.core import has @@ -97,7 +96,7 @@ def eval(cls, expr): elif isinstance(expr, (minus, plus)): return cls(expr, evaluate=False) - elif isinstance(expr, Indexed) and isinstance(expr.base, BasicMapping): + elif isinstance(expr, Indexed) and isinstance(expr.base, SymbolicMapping): return cls(expr, evaluate=False) elif not has(expr, types): if expr.is_number: @@ -106,13 +105,23 @@ def eval(cls, expr): elif isinstance(expr, Expr): x = Symbol(cls.coordinate, real=True) if cls.logical: - M = expr.atoms(Mapping) - if len(M)>0: - M = list(M)[0] - expr_primes = [diff(expr, M[i]) for i in range(M.pdim)] - Jj = Jacobian(M)[:,cls.grad_index] - expr_prime = sum([ei*Jji for ei,Jji in zip(expr_primes, Jj)]) - return expr_prime + diff(expr, x) + # Chain rule here indexes the mapping (`M[i]`), so only an + # indexable point mapping qualifies. A MultiPatchMapping + # (Basic.__new__, no `_shape`) is not indexable and is + # lowered elsewhere. + M = [m for m in expr.atoms(SymbolicMapping) + if getattr(m, '_shape', None) is not None] + if M: + # Chain rule through *every* mapping whose components + # appear in `expr` (D-1). Each contributes + # Sum_i (d expr / d m[i]) * Jacobian(m)[i, k], + # k = cls.grad_index; summation order is irrelevant. + total = diff(expr, x) + for m in M: + Jj = Jacobian(m)[:, cls.grad_index] + total += sum(diff(expr, m[i]) * Jj[i] + for i in range(m.pdim)) + return total return diff(expr, x) @@ -1288,5 +1297,5 @@ def get_max_logical_partial_derivatives(expr, F=None): if v > d[k]: d[k] = v return d -from .mapping import Mapping, Jacobian +from .mapping import Jacobian, SymbolicMapping diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py index 4fe4ebce..d7f14ac7 100644 --- a/sympde/topology/domain.py +++ b/sympde/topology/domain.py @@ -17,16 +17,24 @@ from sympy.tensor import IndexedBase, Indexed from sympy.core import Add, Mul, Pow from sympy.core.expr import AtomicExpr - -from sympde.old_sympy_utilities import is_sequence, with_metaclass -from sympde.core.basic import CalculusFunction +from sympy.utilities.iterables import is_sequence + +from sympde.core.basic import ( + CalculusFunction, + _is_none_argument, + _new_basic, + _restore_container_argument, + _restore_int_argument, + _restore_optional_argument, + _restore_string_argument, +) from .basic import BasicDomain, InteriorDomain, Boundary, Union, Connectivity from .basic import Interval, Interface, CornerBoundary, CornerInterface from .basic import ProductDomain # TODO fix circular dependency between domain and mapping if TYPE_CHECKING: - from sympde.topology.mapping import Mapping + from sympde.topology.mapping import SymbolicMapping # TODO add pdim iterable_types = (tuple, list, Tuple, Union) @@ -41,12 +49,12 @@ class Domain(BasicDomain): name and connectivity need to be passed. """ - def __new__(cls, name : str, *, + def __new__(cls, name : str, *structural_args, interiors : TypeUnion[Iterable[InteriorDomain], InteriorDomain, None] = None, boundaries : TypeUnion[Iterable[Boundary], Boundary, None] = None, dim : Optional[int] = None, connectivity : Optional[Connectivity] = None, - mapping : Optional[Mapping] = None, + mapping : Optional[SymbolicMapping] = None, logical_domain : Optional[Domain] = None): """ Interiors or connectivity must be given. When the mapping is given @@ -69,6 +77,27 @@ def __new__(cls, name : str, *, logical_domain : Domain or None, optional Logical domain that is mapped to the physical domain """ + if structural_args: + if len(structural_args) != 5: + raise TypeError( + 'A reconstructed Domain expects five structural arguments' + ) + if any(value is not None for value in ( + interiors, boundaries, dim, connectivity, + mapping, logical_domain)): + raise TypeError( + 'Structural and keyword Domain arguments cannot be mixed' + ) + interiors, boundaries, mapping, logical_domain, interfaces = structural_args + boundaries = _restore_optional_argument(boundaries) + mapping = _restore_optional_argument(mapping) + logical_domain = _restore_optional_argument(logical_domain) + interfaces = _restore_container_argument(interfaces) + connectivity = Connectivity({str(i.name): i for i in interfaces}) + + name = _restore_string_argument(name) + dim = _restore_int_argument(dim) + # ... if not isinstance(name, str): raise TypeError('> name must be a string') @@ -149,7 +178,8 @@ def __new__(cls, name : str, *, # ... boundaries = Union(*boundaries) - obj = Basic.__new__(cls, name, interiors, boundaries, mapping) + obj = _new_basic(cls, name, interiors, boundaries, mapping, + logical_domain, tuple(connectivity.values())) obj._connectivity = connectivity obj._corners = None obj._dtype = dtype @@ -159,7 +189,7 @@ def __new__(cls, name : str, *, @property def name(self) -> str: - return self.args[0] + return self.args[0].name @property def interior(self) -> TypeUnion[Union, InteriorDomain]: @@ -171,12 +201,14 @@ def interior(self) -> TypeUnion[Union, InteriorDomain]: def boundary(self) -> TypeUnion[Union, Boundary]: """Either a Union object containing the boundaries or just a boundary if there is only one""" - return self.args[2] + boundary = self.args[2] + return None if _is_none_argument(boundary) else boundary @property - def mapping(self) -> Optional[Mapping]: + def mapping(self) -> Optional[SymbolicMapping]: """The mapping that maps the logical domain to the physical domain""" - return self.args[3] + mapping = self.args[3] + return None if _is_none_argument(mapping) else mapping @property def subdomains(self) -> tuple: @@ -375,7 +407,7 @@ def from_file(cls, filename): if ext != '.h5': raise ValueError('> Only h5 files are supported') # ... - from sympde.topology.mapping import Mapping + from sympde.topology.mapping import SymbolicMapping h5 = h5py.File( filename, mode='r' ) yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader ) @@ -398,7 +430,7 @@ def from_file(cls, filename): constructors = [globals()[dt['type']] for dt in dtype] interiors = [cs(i['name'], **dt['parameters']) for cs,i,dt in zip(constructors, d_interior, dtype)] - mappings = [Mapping(I['mapping'], dim=dim) if I.get('mapping', "None") != "None" else None for I in d_interior] + mappings = [SymbolicMapping(I['mapping'], dim=dim) if I.get('mapping', "None") != "None" else None for I in d_interior] domains = [mapping(i) if mapping else i for i,mapping in zip(interiors, mappings)] patch_index = {I.name:ind for ind,I in enumerate(interiors)} @@ -810,8 +842,21 @@ class NCubeInterior(InteriorDomain): def __new__(cls, name, dim=None, dtype=None, min_coords=None, max_coords=None, mapping=None, logical_domain=None): - obj = InteriorDomain.__new__(cls, name, dim=dim, dtype=dtype, - mapping=mapping, logical_domain=logical_domain) + name = _restore_string_argument(name) + dim = _restore_int_argument(dim) + dtype = _restore_container_argument(dtype) + min_coords = _restore_container_argument(min_coords) + max_coords = _restore_container_argument(max_coords) + mapping = _restore_optional_argument(mapping) + logical_domain = _restore_optional_argument(logical_domain) + + obj = _new_basic(cls, name, dim, dtype, min_coords, max_coords, + mapping, logical_domain) + obj._dim = dim + obj._target = None + obj._dtype = dtype + obj._mapping = mapping + obj._logical_domain = logical_domain obj._min_coords = min_coords obj._max_coords = max_coords @@ -838,6 +883,14 @@ def min_coords(self): def max_coords(self): return self._max_coords + @property + def mapping(self): + return self._mapping + + @property + def logical_domain(self): + return self._logical_domain + @property def boundary(self): return self._boundary @@ -865,7 +918,17 @@ def get_boundary(self, axis=None, ext=None): # class NCube(Domain): - def __new__(cls, name, dim, min_coords, max_coords): + def __new__(cls, name, dim, min_coords, max_coords, *structural_args): + + if structural_args: + return Domain.__new__( + cls, name, dim, min_coords, max_coords, *structural_args + ) + + name = _restore_string_argument(name) + dim = _restore_int_argument(dim) + min_coords = _restore_container_argument(min_coords) + max_coords = _restore_container_argument(max_coords) assert isinstance(name, str) assert isinstance(dim, (int, Integer)) @@ -950,7 +1013,9 @@ def max_coords(self): #============================================================================== class Line(NCube): - def __new__(cls, name='Line', bounds=(0, 1)): + def __new__(cls, name='Line', bounds=(0, 1), *structural_args): + if structural_args: + return Domain.__new__(cls, name, bounds, *structural_args) dim = 1 min_coords = (bounds[0],) max_coords = (bounds[1],) @@ -963,7 +1028,12 @@ def bounds(self): #============================================================================== class Square(NCube): - def __new__(cls, name='Square', bounds1=(0, 1), bounds2=(0, 1)): + def __new__(cls, name='Square', bounds1=(0, 1), bounds2=(0, 1), + *structural_args): + if structural_args: + return Domain.__new__( + cls, name, bounds1, bounds2, *structural_args + ) dim = 2 min_coords = (bounds1[0], bounds2[0]) max_coords = (bounds1[1], bounds2[1]) @@ -980,7 +1050,12 @@ def bounds2(self): #============================================================================== class Cube(NCube): - def __new__(cls, name='Cube', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1)): + def __new__(cls, name='Cube', bounds1=(0, 1), bounds2=(0, 1), + bounds3=(0, 1), *structural_args): + if structural_args: + return Domain.__new__( + cls, name, bounds1, bounds2, bounds3, *structural_args + ) dim = 3 min_coords = (bounds1[0], bounds2[0], bounds3[0]) max_coords = (bounds1[1], bounds2[1], bounds3[1]) @@ -1015,7 +1090,7 @@ class TangentVector(BoundaryVector): pass #============================================================================== -class ElementDomain(with_metaclass(Singleton, Basic)): +class ElementDomain(Basic, metaclass=Singleton): pass #============================================================================== diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 585dd85e..eabf1045 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -1,5 +1,6 @@ # coding: utf-8 -from abc import ABC, abstractmethod +import warnings +from abc import ABC, ABCMeta, abstractmethod from sympy import Indexed, IndexedBase, Idx from sympy import Matrix, ImmutableDenseMatrix from sympy import Function, Expr @@ -9,16 +10,17 @@ from sympy.core import Symbol,Integer from sympy.core import Add, Mul, Pow from sympy.core.numbers import ImaginaryUnit -from sympy.core.containers import Tuple +from sympy.core.containers import Dict, Tuple from sympy import S from sympy import sqrt, symbols from sympy.core.exprtools import factor_terms from sympy.polys.polytools import parallel_poly_from_expr +from sympy.utilities.iterables import is_sequence from sympde.core import Constant -from sympde.core.basic import BasicMapping from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery +from sympde.core.basic import _new_basic 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 @@ -44,9 +46,12 @@ # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping __all__ = ( + 'AnalyticMapping', 'BasicCallableMapping', 'Contravariant', 'Covariant', + 'DefinedMapping', + 'DiscreteMapping', 'InterfaceMapping', 'InverseMapping', 'Jacobian', @@ -58,7 +63,9 @@ 'MappingApplication', 'MultiPatchMapping', 'PullBack', + 'StructuralMapping', 'SymbolicExpr', + 'SymbolicMapping', 'SymbolicWeightedVolume', 'get_logical_test_function', ) @@ -95,12 +102,19 @@ def numpy_to_native_python(a): #============================================================================== class BasicCallableMapping(ABC): """ - Transformation of coordinates, which can be evaluated. + Interface for evaluating a coordinate transformation at a point. F: R^l -> R^p F(eta) = x with l <= p + + This has no symbolic identity of its own (no name, not usable in a domain + or form) -- it is the pure evaluation contract that a numeric backend like + psydac's ``SplineCallableMapping`` implements directly. Contrast with + ``SymbolicMapping.__call__``, which maps a *domain* to build a symbolic + ``MappedDomain`` rather than evaluating a point; ``DefinedMapping`` + combines both. """ @abstractmethod def __call__(self, *eta): @@ -139,30 +153,50 @@ def pdim(self): """ #============================================================================== -class Mapping(BasicMapping): +class SymbolicMapping(IndexedBase): """ - Represents a Mapping object. - - Examples - + Common root of the unified mapping hierarchy: an undefined symbolic + transformation of coordinates, identified by a name and a pair of + dimensions (logical ``ldim`` to physical ``pdim``). Callable on a *domain* + (-> a symbolic ``MappedDomain``) but not point-evaluable -- that is + ``DefinedMapping``'s job. + + Construct one directly for a mapping you only need symbolically: + ``SymbolicMapping('F', dim=2)``. Analytic mappings (carrying ``_expressions``) + are ``AnalyticMapping`` subclasses; the ``Mapping`` name (WP06d-4c) is a + deprecated shell over this class. """ - _expressions = None # used for analytical mapping + + _expressions = None _jac = None _inv_jac = None _constants = None + _metric = None + _metric_det = None _callable_map = None _ldim = None _pdim = None + _is_minus = None + _is_plus = None def __new__(cls, name, dim=None, **kwargs): + # `SymbolicMapping.__new__` runs before `object.__new__`'s abstract-class + # guard would fire, so replicate it here: `DefinedMapping('F')` / + # `StructuralMapping('F')` must still raise `TypeError` mentioning the + # unimplemented abstract methods, not trip an assertion below. + if getattr(cls, '__abstractmethods__', frozenset()): + raise TypeError( + "Can't instantiate abstract class {} with abstract methods {}" + .format(cls.__name__, ', '.join(sorted(cls.__abstractmethods__)))) + ldim = kwargs.pop('ldim', cls._ldim) pdim = kwargs.pop('pdim', cls._pdim) coordinates = kwargs.pop('coordinates', None) evaluate = kwargs.pop('evaluate', True) dims = [dim, ldim, pdim] - for i,d in enumerate(dims): + for i, d in enumerate(dims): if isinstance(d, (tuple, list, Tuple, Matrix, ImmutableDenseMatrix)): if not len(d) == 1: raise ValueError('> Expecting a tuple, list, Tuple of length 1') @@ -178,227 +212,235 @@ def __new__(cls, name, dim=None, **kwargs): ldim = dim pdim = dim - obj = IndexedBase.__new__(cls, name, shape=pdim) if not evaluate: return obj if coordinates is None: - _coordinates = [Symbol(name, real=True) for name in ['x', 'y', 'z'][:pdim]] + _coordinates = [Symbol(c, real=True) for c in ['x', 'y', 'z'][:pdim]] else: if not isinstance(coordinates, (list, tuple, Tuple)): raise TypeError('> Expecting list, tuple, Tuple') - for a in coordinates: if not isinstance(a, (str, Symbol)): raise TypeError('> Expecting str or Symbol') - _coordinates = [Symbol(u, real=True) for u in coordinates] obj._name = name obj._ldim = ldim obj._pdim = pdim obj._coordinates = tuple(_coordinates) + # `jacobian=` lets a caller inject a pre-built symbolic Jacobian + # instead of the default JacobianSymbol(obj); preserved from the old + # Mapping.__new__ constructor contract. obj._jacobian = kwargs.pop('jacobian', JacobianSymbol(obj)) obj._is_minus = None obj._is_plus = None - lcoords_names = ['x1', 'x2', 'x3'][:ldim] + 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_names = ['x', 'y', 'z'][:pdim] - - # ... - args = [] - for i in coords_names: - x = obj._expressions[i] - x = sympify(x) - args.append(x) - - args = Tuple(*args) - # ... - zero_coords = ['x1', 'x2', 'x3'][ldim:] - - for i in zero_coords: - x = sympify(i) - args = args.subs(x,0) - # ... - # 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)) - - args = [obj[i] for i in range(pdim)] - exprs = obj._expressions - subs = list(zip(_coordinates, exprs)) - - if obj._jac is None and obj._inv_jac is None: - obj._jac = Jacobian(obj).subs(list(zip(args, exprs))) - obj._inv_jac = obj._jac.inv() if pdim == ldim else None - elif obj._inv_jac is None: - obj._jac = ImmutableDenseMatrix(sympify(obj._jac)).subs(subs) - obj._inv_jac = obj._jac.inv() if pdim == ldim else None - - elif obj._jac is None: - obj._inv_jac = ImmutableDenseMatrix(sympify(obj._inv_jac)).subs(subs) - obj._jac = obj._inv_jac.inv() - else: - obj._jac = ImmutableDenseMatrix(sympify(obj._jac)).subs(subs) - obj._inv_jac = ImmutableDenseMatrix(sympify(obj._inv_jac)).subs(subs) - - else: - obj._jac = Jacobian(obj) + obj._logical_coordinates = Tuple(*lcoords_symbols_real) - obj._metric = obj._jac.T*obj._jac - obj._metric_det = obj._metric.det() + # Undefined mapping: symbolic Jacobian / metric straight from + # Jacobian(obj). Mapping.__new__ overrides this for the subclasses that + # carry `_expressions`. + if cls._expressions is None: + obj._jac = Jacobian(obj) + obj._metric = obj._jac.T * obj._jac + obj._metric_det = obj._metric.det() return obj + # Applying the mapping to a logical domain returns a mapped domain. + def __call__(self, domain): + assert isinstance(domain, BasicDomain) + assert domain.logical_domain is None + assert domain.dim == self.ldim + return MappedDomain(self, domain) + #-------------------------------------------------------------------------- - # Callable mapping + # Callable mapping (a discrete point-evaluable mapping attached explicitly) #-------------------------------------------------------------------------- def get_callable_mapping(self): - if self._callable_map is None: - if self._expressions is None: - msg = 'Cannot generate callable mapping without analytical expressions. '\ - 'A user-defined callable mapping of type `BasicCallableMapping` '\ - 'can be provided using the method `set_callable_mapping`.' - raise ValueError(msg) - - from sympde.topology.callable_mapping import CallableMapping - self._callable_map = CallableMapping(self) - - return self._callable_map + # An explicitly attached callable (set_callable_mapping) always wins. + if self._callable_map is not None: + return self._callable_map + + if self._expressions is None: + raise ValueError( + 'Cannot generate a callable mapping without analytical ' + 'expressions. Give this mapping a point-evaluable identity ' + 'up front instead, e.g. `DiscreteMapping(F, name)` (or ' + '`F.to_defined_mapping(name)` if `F` is a psydac ' + '`SplineCallableMapping`).') + + # Reachable only for a bare subclass that carries `_expressions` but is + # not an `AnalyticMapping` (whose override returns `self`). Since WP06c + # such a class is a mistake -- there is no `CallableMapping` to build. + raise TypeError( + f'{type(self).__name__} carries analytical expressions but is not ' + 'an AnalyticMapping. Subclass AnalyticMapping for a point-evaluable ' + 'analytic mapping; an AnalyticMapping instance is its own callable ' + 'mapping.') def set_callable_mapping(self, F): + """ + Attach a point-evaluable callable to this undefined symbolic mapping. + + .. deprecated:: 0.18 + Mutating an already-constructed symbolic mapping is unsafe: + sympy interns ``SymbolicMapping`` instances, so every other holder + of a mapping with the same name and dimensions observes the + change. Build the identity up front instead -- + ``DiscreteMapping(F, name)`` -- and map the logical domain with + that. The method still works and will be removed in a later step. + Parameters + ---------- + F : BasicCallableMapping + The callable to attach. + + Raises + ------ + TypeError + If `F` is not a `BasicCallableMapping`. (Also raised + unconditionally by `AnalyticMapping`, `DiscreteMapping` and + `StructuralMapping`, which override this method.) + + Note the deprecation warning above is emitted *before* this + check, deliberately, so that a caller on the deprecated path + learns it is deprecated even when their argument is also wrong. + Under warnings-as-errors (``-W error::DeprecationWarning``) that + means a bad `F` surfaces as the `DeprecationWarning`, not as this + `TypeError`. + + Examples + -------- + >>> F = SymbolicMapping('F', dim=2) + >>> F.set_callable_mapping(F_h) # deprecated + >>> F.get_callable_mapping() is F_h + True + + The replacement, which builds the identity up front: + + >>> G = DiscreteMapping(F_h, 'G') + >>> G.get_callable_mapping() is F_h + True + """ + warnings.warn( + 'SymbolicMapping.set_callable_mapping is deprecated: it mutates a ' + 'sympy-interned object, so other holders of the same symbolic ' + 'mapping silently become point-evaluable too. Use ' + 'DiscreteMapping(F, name) instead (or F.to_defined_mapping(name) ' + 'if F is a psydac SplineCallableMapping).', + DeprecationWarning, stacklevel=2) if not isinstance(F, BasicCallableMapping): raise TypeError( f'F must be a BasicCallableMapping, got {type(F)} instead') - self._callable_map = F #-------------------------------------------------------------------------- @property - def name( self ): + def name(self): return self._name @property - def ldim( self ): + def ldim(self): return self._ldim @property - def pdim( self ): + def pdim(self): return self._pdim @property - def coordinates( self ): - if self.pdim == 1: - return self._coordinates[0] - else: - return self._coordinates + def coordinates(self): + return self._coordinates[0] if self.pdim == 1 else self._coordinates @property - def logical_coordinates( self ): - if self.ldim == 1: - return self._logical_coordinates[0] - else: - return self._logical_coordinates - - # Applying the mapping to a logical domain returns a mapped domain - def __call__(self, domain): - assert isinstance(domain, BasicDomain) - assert domain.logical_domain is None - assert domain.dim == self.ldim - return MappedDomain(self, domain) + def logical_coordinates(self): + return self._logical_coordinates[0] if self.ldim == 1 else self._logical_coordinates @property - def jacobian( self ): + def jacobian_symbol(self): + """ Symbolic JacobianSymbol of this mapping. """ return self._jacobian @property - def det_jacobian( self ): - return self.jacobian.det() + def det_jacobian(self): + return self.jacobian_symbol.det() @property - def is_analytical( self ): - return not( self._expressions is None ) + def is_analytical(self): + return self._expressions is not None @property - def expressions( self ): + def expressions(self): return self._expressions @property - def jacobian_expr( self ): + def constants(self): + return self._constants + + @property + def jacobian_expr(self): return self._jac @property - def jacobian_inv_expr( self ): + def jacobian_inv_expr(self): if not self.is_analytical and self._inv_jac is None: self._inv_jac = self.jacobian_expr.inv() return self._inv_jac @property - def metric_expr( self ): + def metric_expr(self): return self._metric @property - def metric_det_expr( self ): + def metric_det_expr(self): return self._metric_det @property - def constants( self ): - return self._constants - - @property - def is_minus( self ): + def is_minus(self): return self._is_minus @property - def is_plus( self ): + def is_plus(self): return self._is_plus - def set_plus_minus( self, **kwargs): + def set_plus_minus(self, **kwargs): minus = kwargs.pop('minus', False) plus = kwargs.pop('plus', False) assert plus is not minus - self._is_plus = plus self._is_minus = minus def copy(self): - obj = Mapping(self.name, - ldim=self.ldim, - pdim=self.pdim, - evaluate=False) - + # type(self), not SymbolicMapping, so a copy keeps its concrete class + # (and its point-evaluation capability, for AnalyticMapping). Safe + # because evaluate=False short-circuits __new__ right after + # IndexedBase.__new__, before the `_expressions` branch. + # + # Suppress the bare-`Mapping` deprecation warning: this is an internal + # reconstruction, not a new user construction (InterfaceMapping.__new__ + # copies its legs). + with warnings.catch_warnings(): + warnings.simplefilter('ignore', DeprecationWarning) + obj = type(self)(self.name, ldim=self.ldim, pdim=self.pdim, + evaluate=False) obj._name = self.name obj._ldim = self.ldim obj._pdim = self.pdim - obj._coordinates = self.coordinates + obj._coordinates = self._coordinates obj._jacobian = JacobianSymbol(obj) - obj._logical_coordinates = self.logical_coordinates + obj._logical_coordinates = self._logical_coordinates obj._expressions = self._expressions obj._constants = self._constants obj._jac = self._jac obj._inv_jac = self._inv_jac obj._metric = self._metric obj._metric_det = self._metric_det - obj.__callable_map = self._callable_map + obj._callable_map = self._callable_map obj._is_plus = self._is_plus obj._is_minus = self._is_minus return obj @@ -412,25 +454,753 @@ def _eval_subs(self, old, new): return self def _sympystr(self, printer): - sstr = printer.doprint - return sstr(self.name) + return printer.doprint(self.name) + +#============================================================================== +class _MappingABCMeta(ABCMeta, type(SymbolicMapping)): + """ + Metaclass merging ``abc.ABCMeta`` with sympy's metaclass + (``ManagedProperties`` in sympy 1.9). + + A mapping class derived from ``IndexedBase`` (via ``SymbolicMapping``) + already carries sympy's metaclass; declaring ``@abstractmethod`` members on + it additionally requires ``ABCMeta``. Python rejects a class whose metaclass + is not a subclass of every base's metaclass, so the two are merged here once + and reused by ``DefinedMapping``, ``StructuralMapping``, and + ``AnalyticMapping``. ``type(...)`` is used instead of importing the name so + this keeps working if a future sympy renames its metaclass. + """ + +#============================================================================== +class DefinedMapping(SymbolicMapping, BasicCallableMapping, metaclass=_MappingABCMeta): + """ + Abstract base class for *point-evaluable* mappings. + + F: R^l -> R^p , F(eta) = x , with l <= p + + The two concrete subclasses are ``AnalyticMapping`` (lambdifies its own + ``_expressions``) and ``DiscreteMapping`` (delegates to a wrapped external + :class:`BasicCallableMapping`, e.g. psydac's ``SplineCallableMapping``). A + bare :class:`BasicCallableMapping` has no symbolic identity of its own -- + it must be wrapped (``DiscreteMapping(F, name)`` / + ``F.to_defined_mapping(name)``) to act as a ``DefinedMapping``. The + point-evaluation interface (``__call__``, ``jacobian``, ``jacobian_inv``, + ``metric``, ``metric_det``, ``ldim``, ``pdim``) is inherited verbatim from + ``BasicCallableMapping``; + every one must be implemented for a subclass to be instantiable, which is + what makes the two subclasses interchangeable wherever a + :class:`BasicCallableMapping` is expected. + + ``DefinedMapping`` itself is abstract, but as a convenience + ``DefinedMapping(callable_mapping, name, dim=...)`` -- a first positional + :class:`BasicCallableMapping` -- builds a :class:`DiscreteMapping` wrapping + that callable. + """ + + def __new__(cls, *args, **kwargs): + if (cls is DefinedMapping and args + and isinstance(args[0], BasicCallableMapping)): + return DiscreteMapping(*args, **kwargs) + return super().__new__(cls, *args, **kwargs) + +#============================================================================== +class StructuralMapping(SymbolicMapping, metaclass=_MappingABCMeta): + """ + Abstract base class for *symbolic, non-point-evaluable* mappings. + + These objects (``InverseMapping``, ``InterfaceMapping``, + ``MultiPatchMapping``) describe how patches are related or assembled and + only make sense symbolically: they stay callable on a domain (returning a + symbolic mapped domain, like any ``SymbolicMapping``) but reject point + evaluation. + + Since WP06d-4a a ``StructuralMapping`` no longer inherits ``Mapping`` -- it + owns the small symbolic surface it needs (``name``, ``is_plus``/``is_minus``, + coordinates, ``__call__``, hashing) directly here. Each concrete subclass + sets ``_name`` / ``_ldim`` / ``_pdim`` (and ``_coordinates`` / + ``_logical_coordinates`` where meaningful) in its own ``__new__`` and + implements ``ldim`` / ``pdim``. + """ + + _is_minus = None + _is_plus = None + + def __call__(self, *args, domain=None): + """ + Call this structural mapping on a domain (positional or as the + ``domain`` keyword, matching ``SymbolicMapping``'s signature) to get a + symbolic mapped domain. + + Parameters + ---------- + domain : BasicDomain + The logical domain to map. + + Returns + ------- + MappedDomain + + Raises + ------ + TypeError + If called with anything other than exactly one ``BasicDomain`` + argument -- a ``StructuralMapping`` is symbolic and not + point-evaluable. + """ + # Same domain-vs-point dispatch as AnalyticMapping.__call__; the + # terminal action here is the domain call (Mapping.__call__'s body, + # reproduced -- Mapping is no longer in the MRO to super() into). + if domain is None and len(args) == 1 and isinstance(args[0], BasicDomain): + domain = args[0] + if domain is not None: + assert domain.logical_domain is None + assert domain.dim == self.ldim + return MappedDomain(self, domain) + raise TypeError( + f"{type(self).__name__} is a StructuralMapping: it is symbolic " + "and not point-evaluable.") + + def set_callable_mapping(self, F): + # WP09: a StructuralMapping is symbolic and not point-evaluable (see + # __call__) -- attaching a callable would silently make + # get_callable_mapping() return one, contradicting that. Mirrors + # DiscreteMapping's guard (WP07d). + raise TypeError( + f"{type(self).__name__} is a StructuralMapping: it is symbolic " + "and not point-evaluable, and cannot be given an attached " + "callable.") + + @property + def name(self): + return self._name + + @property + def is_minus(self): + return self._is_minus + + @property + def is_plus(self): + return self._is_plus + + @property + def coordinates(self): + c = self._coordinates + return c[0] if self.pdim == 1 else c + + @property + def logical_coordinates(self): + c = self._logical_coordinates + return c[0] if self.ldim == 1 else c + + @property + def jacobian_symbol(self): + # lazy: multipatch pull-back (PullBack.__new__) asks for this on an + # InterfaceMapping. InverseMapping overrides with a pre-inverted one. + j = getattr(self, '_jacobian', None) + if j is None: + j = self._jacobian = JacobianSymbol(self) + return j + + # A structural mapping carries no analytic `_expressions`; its symbolic + # jacobian / metric are the same `Jacobian(self)` forms that Mapping.__new__ + # produced for it before WP06d-4a severed the base. Each concrete __new__ + # calls _init_symbolic_jacobian(obj) once (straight-line, like the old + # Mapping.__new__ `else` branch -- a lazy property would recurse through + # Jacobian.eval). psydac codegen (api/ast/fem.py) reads `jacobian_expr` on + # multipatch mappings. + _jac = None + _inv_jac = None + _metric = None + _metric_det = None + + @staticmethod + def _init_symbolic_jacobian(obj): + obj._jac = Jacobian(obj) + obj._metric = obj._jac.T * obj._jac + obj._metric_det = obj._metric.det() + + @property + def jacobian_expr(self): + return self._jac + + @property + def jacobian_inv_expr(self): + return self._inv_jac + + @property + def metric_expr(self): + return self._metric + + @property + def metric_det_expr(self): + return self._metric_det + + def _hashable_content(self): + return (type(self).__name__, self._name, self.ldim, self.pdim, + self._coordinates, self._logical_coordinates) + + @property + @abstractmethod + def ldim(self): + """ Number of logical/parametric dimensions. """ + + @property + @abstractmethod + def pdim(self): + """ Number of physical dimensions. """ + +#============================================================================== +def _init_analytic_expressions(obj, **kwargs): + """ + Expand ``obj._expressions`` (set by the class body of an analytic + mapping subclass) into ``_jac`` / ``_inv_jac`` / ``_metric`` / + ``_metric_det`` / ``_constants``, and return ``obj``. ``kwargs`` are the + caller's numeric values for the symbolic constants found in + ``_expressions`` (anything not given numerically becomes a + :class:`~sympde.core.basic.Constant`). + + Free-standing (not a method) because it is shared by two callers: + ``AnalyticMapping.__new__`` (the supported path) and the deprecated + ``Mapping.__new__`` (kept working for legacy subclasses such as + ``class SquareTorus(Mapping)`` in downstream code). Keeping it + free-standing also means ``Mapping`` can be deleted later, by the + remove-aliases work-package, without moving this code again. + """ + # -- analytic subclass: expand `_expressions` into _jac / _inv_jac / + # _metric (SymbolicMapping.__new__ left these to us). -- + ldim, pdim = obj._ldim, obj._pdim + _coordinates = list(obj._coordinates) + lcoords_names = ['x1', 'x2', 'x3'][:ldim] + lcoords_symbols_real = list(obj._logical_coordinates) + lcoords_symbols_real_dict = {n: s for n, s in zip(lcoords_names, lcoords_symbols_real)} + lcoords_general_symbols = [Symbol(i) for i in lcoords_names] + + coords_names = ['x', 'y', 'z'][:pdim] + args = Tuple(*(sympify(obj._expressions[i]) for i in coords_names)) + for i in ['x1', 'x2', 'x3'][ldim:]: + args = args.subs(sympify(i), 0) + + constants = list(set(args.free_symbols) - set(lcoords_general_symbols)) + constants_values = {a.name: Constant(a.name) for a in constants} + constants_values.update(kwargs) + d = {a: numpy_to_native_python(constants_values[a.name]) for a in constants} + args = args.subs(d) + 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)) + + idx = [obj[i] for i in range(pdim)] + exprs = obj._expressions + subs = list(zip(_coordinates, exprs)) + + if obj._jac is None and obj._inv_jac is None: + obj._jac = Jacobian(obj).subs(list(zip(idx, exprs))) + obj._inv_jac = obj._jac.inv() if pdim == ldim else None + elif obj._inv_jac is None: + obj._jac = ImmutableDenseMatrix(sympify(obj._jac)).subs(subs) + obj._inv_jac = obj._jac.inv() if pdim == ldim else None + elif obj._jac is None: + obj._inv_jac = ImmutableDenseMatrix(sympify(obj._inv_jac)).subs(subs) + obj._jac = obj._inv_jac.inv() + else: + obj._jac = ImmutableDenseMatrix(sympify(obj._jac)).subs(subs) + obj._inv_jac = ImmutableDenseMatrix(sympify(obj._inv_jac)).subs(subs) + + obj._metric = obj._jac.T * obj._jac + obj._metric_det = obj._metric.det() + return obj #============================================================================== -class InverseMapping(Mapping): +class Mapping(SymbolicMapping): + """ + Deprecated shell over :class:`SymbolicMapping` (WP06d-4c): constructing a + bare ``Mapping`` emits a ``DeprecationWarning`` -- use ``SymbolicMapping`` + for an undefined mapping, or an ``AnalyticMapping`` subclass for an analytic + one. It is no longer in ``AnalyticMapping``'s MRO (WP D8): the + ``_expressions`` -> symbolic Jacobian / metric machinery now lives in the + free-standing :func:`_init_analytic_expressions`, which ``Mapping.__new__`` + still calls so that legacy subclasses (``class X(Mapping)`` with + ``_expressions``) keep working. ``Mapping`` survives only as a deprecated + name and as that legacy analytic base, until the remove-aliases + work-package deletes it. + + .. deprecated:: WP D8 + Use :class:`SymbolicMapping` (undefined mapping) or + :class:`AnalyticMapping` (analytic mapping) instead. + """ + + def __new__(cls, name, dim=None, **kwargs): + # Warn only on explicit user construction (`name` is a plain ``str``). + # sympy's internal ``func(*args)`` rebuild -- pickling, ``Basic.rebuild``, + # cse/canonicalisation -- re-invokes ``__new__`` with ``name`` already a + # ``Symbol``; that path must stay quiet so a downstream running + # ``filterwarnings=error`` is not tripped by warning spam it did not + # cause. + if cls is Mapping and isinstance(name, str): + warnings.warn( + 'sympde.topology.Mapping is deprecated; use SymbolicMapping ' + '(or an AnalyticMapping subclass for an analytic mapping).', + DeprecationWarning, stacklevel=2) + + evaluate = kwargs.get('evaluate', True) + # `jacobian` goes to SymbolicMapping.__new__ too -- forwarding it here + # both honours it and keeps it out of the `_expressions` constants dict. + prefix_kw = {k: kwargs.pop(k) + for k in ('ldim', 'pdim', 'coordinates', 'evaluate', 'jacobian') + if k in kwargs} + obj = super().__new__(cls, name, dim, **prefix_kw) + + if not evaluate or cls._expressions is None: + return obj + + return _init_analytic_expressions(obj, **kwargs) + + @property + def jacobian(self): + # WP02a froze this name so WP02b/06c could make it the *numeric* + # point-evaluation method on DefinedMapping. The symbolic JacobianSymbol + # lives on `jacobian_symbol`. + warnings.warn( + "Mapping.jacobian (symbolic) is deprecated; use " + "Mapping.jacobian_symbol.", + DeprecationWarning, stacklevel=2) + return self.jacobian_symbol + +#============================================================================== +class AnalyticMapping(DefinedMapping, metaclass=_MappingABCMeta): + """ + Analytic mapping: a concrete :class:`DefinedMapping` whose ``_expressions`` + (set by the subclass body) are expanded into symbolic Jacobian / metric on + construction, and *directly point-evaluable* through the + ``DefinedMapping`` interface. + + Point evaluation is done by the mapping itself: on first use it lambdifies + its own analytic expressions into numpy callables (cached, one quantity at a + time), so an ``AnalyticMapping`` and a psydac ``SplineCallableMapping`` are + interchangeable wherever the point-evaluation interface is expected, and + ``get_callable_mapping()`` returns ``self``. + + A subclass **must** define a non-``None`` ``_expressions`` dict in its class + body (see the gallery in ``analytical_mapping.py``); it is what makes the + mapping point-evaluable. Instantiating a subclass that leaves + ``_expressions`` unset raises ``TypeError`` -- D11: without expressions + there is nothing to lambdify, so the alternative would be a + ``DefinedMapping`` in name only, silently unusable until some unrelated + call fails far from where it was built. Use ``SymbolicMapping(name, + dim=...)`` for an undefined symbolic mapping instead. Setting + ``cls._expressions`` before calling ``super().__new__`` (e.g. from a + factory) still works, since the guard reads ``cls._expressions``. + + Examples + -------- + >>> from sympde.topology import IdentityMapping + >>> F = IdentityMapping('F', dim=2) + >>> F(0.3, 0.4) # point evaluation -> physical coords + (0.3, 0.4) + >>> F.jacobian_symbol # symbolic JacobianSymbol (unchanged) + Jacobian(F) + + Raises + ------ + TypeError + If the class being instantiated -- ``AnalyticMapping`` itself or any + subclass -- has no ``_expressions``, either its own or inherited. + Only *instantiating* such a class raises; *defining* one (e.g. an + intermediate base) is still allowed. + """ + + def __new__(cls, name, dim=None, **kwargs): + # D11: without `_expressions` there is nothing to point-evaluate, so the + # result would be a DefinedMapping in name only (it used to construct and + # then fail far away -- e.g. discretize(Derham) succeeded, assembly died + # in generated code). Checked on `cls`, at instantiation (not in + # __init_subclass__), so an intermediate base without `_expressions` can + # still be *defined* and subclassed. Unconditional on `evaluate`: the only + # evaluate=False caller, SymbolicMapping.copy(), copies an existing + # instance, whose class already passed this check. Placed before + # `super().__new__` so no half-built instance is ever interned by sympy. + if cls._expressions is None: + raise TypeError( + f"{cls.__name__} has no analytical expressions (`_expressions` " + "is None), so it cannot be point-evaluated. Use " + "SymbolicMapping(name, dim=...) for an undefined symbolic " + "mapping, subclass AnalyticMapping with an `_expressions` dict " + "for an analytic one, or DiscreteMapping(F, name) to wrap an " + "external point-evaluable callable.") + evaluate = kwargs.get('evaluate', True) + # `jacobian` goes to SymbolicMapping.__new__ too -- forwarding it here + # both honours it and keeps it out of the `_expressions` constants dict. + prefix_kw = {k: kwargs.pop(k) + for k in ('ldim', 'pdim', 'coordinates', 'evaluate', 'jacobian') + if k in kwargs} + obj = super().__new__(cls, name, dim, **prefix_kw) + + if not evaluate: + return obj + + return _init_analytic_expressions(obj, **kwargs) + + def _ensure_lambdified(self): + """ Check this mapping can be point-evaluated and return the per-quantity + lambdified-callable cache (created empty on first call; entries are + filled in on demand by ``_lambdify``). """ + cache = getattr(self, '_lambdified', None) + if cache is None: + if self._expressions is None: + raise ValueError('Cannot point-evaluate a mapping without ' + 'analytical expressions.') + if self._constants: + raise ValueError( + f'{self.name} has unresolved symbolic constants ' + f'{self._constants}; give them numeric values at ' + 'construction to point-evaluate it.') + self._lambdified = cache = {} + return cache + + def _lambdify(self, key): + """ Lambdify one quantity -- ``'call'``, ``'jacobian'``, + ``'jacobian_inv'``, ``'metric'`` or ``'metric_det'`` -- into a numpy + callable, caching it on the instance. Done one quantity at a time so + that e.g. ``F(eta)`` never builds ``jacobian_inv`` (which needs a + square Jacobian -- surface mappings have ``ldim != pdim``). """ + cache = self._ensure_lambdified() + if key not in cache: + # Lazy import: sympde.utilities.utils imports sympde.topology, + # so a module-level import here would be circular. + from sympde.utilities.utils import lambdify_sympde + v = self.logical_coordinates + if key == 'call': + cache[key] = tuple(lambdify_sympde(v, e) for e in self.expressions) + else: + expr = {'jacobian': self.jacobian_expr, + 'jacobian_inv': self.jacobian_inv_expr, + 'metric': self.metric_expr, + 'metric_det': self.metric_det_expr}[key] + cache[key] = lambdify_sympde(v, expr) + return cache[key] + + def _delegate_point_eval(self, name, *eta): + """ Point-evaluate quantity ``name`` (``'call'`` | ``'jacobian'`` | + ``'jacobian_inv'`` | ``'metric'`` | ``'metric_det'``) at ``eta``, using + this mapping's own lambdified analytic expressions. """ + if name == 'call': + return tuple(f(*eta) for f in self._lambdify('call')) + return self._lambdify(name)(*eta) + + def __call__(self, *args, domain=None): + # A single BasicDomain (positional or `domain=`) -> symbolic + # MappedDomain (unchanged Mapping behaviour). Anything else -> point + # evaluation on logical coordinates. An unexpected keyword raises a + # natural TypeError. Same dispatch shape as StructuralMapping.__call__. + if domain is None and len(args) == 1 and isinstance(args[0], BasicDomain): + domain, args = args[0], () + if domain is not None: + return super().__call__(domain) + return self._delegate_point_eval('call', *args) + + def jacobian(self, *eta): + """ Jacobian matrix evaluated at the logical point(s) ``eta``. """ + return self._delegate_point_eval('jacobian', *eta) + + def jacobian_inv(self, *eta): + """ Inverse Jacobian matrix evaluated at the logical point(s) ``eta``. """ + return self._delegate_point_eval('jacobian_inv', *eta) + + def metric(self, *eta): + """ Metric tensor evaluated at the logical point(s) ``eta``. """ + return self._delegate_point_eval('metric', *eta) + + def metric_det(self, *eta): + """ Determinant of the metric tensor at the logical point(s) ``eta``. """ + return self._delegate_point_eval('metric_det', *eta) + + def get_callable_mapping(self): + # An AnalyticMapping *is* its own callable mapping -- guaranteed since + # D11: construction requires `_expressions`. + return self + + def set_callable_mapping(self, F): + # D2/WP13: an AnalyticMapping is its own callable mapping (self is + # always point-evaluable); attaching an external one would leave + # is_analytical == True while assembly silently kept using the + # analytic formulas instead of the attached callable. Mirrors + # DiscreteMapping's guard (WP07d) and StructuralMapping's (WP09). + # Use DiscreteMapping(F, name) instead (or, if F is a psydac + # SplineCallableMapping specifically, its convenience shorthand + # F.to_defined_mapping(name) -- not available on a generic + # BasicCallableMapping). + raise TypeError( + f"{type(self).__name__} is an AnalyticMapping: it is already its " + "own callable mapping, and cannot be given a different attached " + "callable. Wrap the callable instead: DiscreteMapping(F, name) " + "(or F.to_defined_mapping(name) if F is a psydac " + "SplineCallableMapping).") + + # ldim / pdim: the concrete properties inherited from SymbolicMapping + # satisfy the DefinedMapping / BasicCallableMapping abstract members. + +#============================================================================== +class DiscreteMapping(DefinedMapping, metaclass=_MappingABCMeta): + """ + A concrete :class:`DefinedMapping` whose geometry is supplied by an + *external* point-evaluable mapping (a psydac ``SplineCallableMapping`` / + ``NurbsCallableMapping``, or any :class:`BasicCallableMapping`). + + It plays the same role as :class:`AnalyticMapping` -- a symbolic mapping + that is also point-evaluable -- but where ``AnalyticMapping`` lambdifies its + own ``_expressions``, ``DiscreteMapping`` *delegates* every point call to + the wrapped callable. Use it to give a discrete geometry (e.g. a spline + approximation) a first-class symbolic identity without mutating some other + mapping via ``set_callable_mapping`` (deprecated) -- and unlike that + pattern, the wrapped callable here cannot itself be reattached after + construction (:meth:`set_callable_mapping` raises ``TypeError``; it is + part of :meth:`_hashable_content`):: + + >>> from psydac.mapping.discrete import SplineCallableMapping + >>> F_h = SplineCallableMapping.from_mapping(V, F) # a BasicCallableMapping + >>> G = DiscreteMapping(F_h, name='G', dim=2) # a fresh DefinedMapping + >>> G.get_callable_mapping() is F_h + True + >>> G(domain_log) # callable on a domain -> MappedDomain + >>> G.jacobian(0.3, 0.4) # point call -> delegated to F_h + + ``is_analytical`` is ``False`` (there are no closed-form ``_expressions``), + so psydac's assembly AST takes its grid-evaluation branch for a + ``DiscreteMapping``-carried domain (as it does for a bare + :class:`SymbolicMapping` from ``Domain.from_file``), evaluating the wrapped + callable on the quadrature grid rather than an analytic Jacobian. That + callable must be a ``SplineCallableMapping`` / ``NurbsCallableMapping`` + for assembly today; an end-to-end spline-mapped solve through + ``DiscreteMapping`` is verified in WP07b. Point evaluation, plotting and + symbolic use work with any ``BasicCallableMapping``. + + Two same-named ``DiscreteMapping``s wrapping *different* callables are + distinct (equality / hashing include the wrapped callable), so a per-patch + spline geometry is never silently conflated with another. Still, give each + discrete geometry a distinct ``name`` -- the name is its symbolic identity. + + A ``DiscreteMapping`` has no symbolic sub-structure -- the wrapped callable + is not in ``.args`` and the name / dims do not simplify -- so ``subs`` / + ``xreplace`` and any in-process ``func(*args)`` rebuild (``sympy.simplify``, + ``Basic.rebuild``, cse walking an expression that contains ``self[i]``) are + the identity (see :meth:`func`). ``pickle`` / ``copy.deepcopy`` still cannot + round-trip it -- they reconstruct through ``__new__`` with the symbolic args + only, and the runtime callable is gone; use :meth:`copy` or rebuild from the + original callable. (``InverseMapping`` / ``InterfaceMapping`` / + ``MultiPatchMapping`` share the pickle limitation.) + + Parameters + ---------- + callable_mapping : BasicCallableMapping + The point-evaluable mapping providing the geometry. It is authoritative + for ``ldim`` / ``pdim``. + name : str, optional + Symbolic name. Falls back to ``callable_mapping.name`` when not given; a + ``ValueError`` is raised if neither yields a non-empty name. + ldim, pdim, dim : int, optional + May only *confirm* ``callable_mapping.ldim`` / ``.pdim`` -- a mismatch + raises ``ValueError``. ``dim`` sets both. No other keyword is accepted. + + Raises + ------ + TypeError + If ``callable_mapping`` is not a :class:`BasicCallableMapping`, or an + unexpected keyword is passed. + ValueError + If no non-empty name can be determined, or a supplied ``dim`` / ``ldim`` + / ``pdim`` contradicts the wrapped callable. + """ + + def __new__(cls, callable_mapping, name=None, **kwargs): + if not isinstance(callable_mapping, BasicCallableMapping): + raise TypeError( + 'DiscreteMapping wraps a runtime BasicCallableMapping and ' + 'cannot be reconstructed from symbolic args alone (got ' + f'{type(callable_mapping).__name__}). This happens if sympy ' + 'rebuilds it via func(*args) / pickle / deepcopy; rebuild it ' + 'from the original callable instead. (InverseMapping / ' + 'InterfaceMapping / MultiPatchMapping share this limitation.)') + + if name is None: + name = getattr(callable_mapping, 'name', None) + if not name: + raise ValueError('DiscreteMapping needs a non-empty name (pass ' + 'name=..., or wrap a callable that has a .name)') + + # The wrapped callable is authoritative for the dimensions; a `dim` / + # `ldim` / `pdim` argument may only confirm them, never override. + dim = kwargs.pop('dim', None) + ldim = kwargs.pop('ldim', dim) + pdim = kwargs.pop('pdim', dim) + c_ldim, c_pdim = callable_mapping.ldim, callable_mapping.pdim + if ldim is not None and ldim != c_ldim: + raise ValueError(f'ldim={ldim} conflicts with ' + f'callable_mapping.ldim={c_ldim}') + if pdim is not None and pdim != c_pdim: + raise ValueError(f'pdim={pdim} conflicts with ' + f'callable_mapping.pdim={c_pdim}') + if kwargs: + raise TypeError(f'DiscreteMapping got unexpected keyword(s) ' + f'{list(kwargs)}') + + # Build the IndexedBase directly and attach `_callable_map` BEFORE the + # symbolic Jacobian / metric machinery, so every Indexed(obj, i) it + # bakes in carries `obj`'s final identity -- which includes the callable + # (see `_hashable_content`). Same construction shape as InverseMapping / + # InterfaceMapping; the Jacobian tail matches + # StructuralMapping._init_symbolic_jacobian (inlined -- a DefinedMapping + # calling a StructuralMapping staticmethod would be odd). + obj = IndexedBase.__new__(cls, name, shape=c_pdim) + obj._name = name + obj._ldim = c_ldim + obj._pdim = c_pdim + obj._coordinates = tuple(Symbol(u, real=True) + for u in ['x', 'y', 'z'][:c_pdim]) + obj._logical_coordinates = Tuple(*(Symbol(u, real=True) + for u in ['x1', 'x2', 'x3'][:c_ldim])) + obj._is_minus = None + obj._is_plus = None + obj._callable_map = callable_mapping + obj._jacobian = JacobianSymbol(obj) + obj._jac = Jacobian(obj) + obj._metric = obj._jac.T * obj._jac + obj._metric_det = obj._metric.det() + return obj + + #-------------------------------------------------------------------------- + def __call__(self, *args, domain=None): + # A single BasicDomain (positional or `domain=`) -> symbolic + # MappedDomain; anything else -> point evaluation, delegated to the + # wrapped callable. Same dispatch shape as AnalyticMapping.__call__. + if domain is None and len(args) == 1 and isinstance(args[0], BasicDomain): + domain, args = args[0], () + if domain is not None: + return super().__call__(domain) + return self.get_callable_mapping()(*args) + + def jacobian(self, *eta): + """ Jacobian matrix at the logical point(s) ``eta`` (delegated). """ + return self.get_callable_mapping().jacobian(*eta) + + def jacobian_inv(self, *eta): + """ Inverse Jacobian matrix at the logical point(s) ``eta`` (delegated). """ + return self.get_callable_mapping().jacobian_inv(*eta) + + def metric(self, *eta): + """ Metric tensor at the logical point(s) ``eta`` (delegated). """ + return self.get_callable_mapping().metric(*eta) + + def metric_det(self, *eta): + """ Determinant of the metric tensor at ``eta`` (delegated). """ + return self.get_callable_mapping().metric_det(*eta) + + def get_callable_mapping(self): + # Unlike AnalyticMapping, a DiscreteMapping is never its own callable. + if self._callable_map is None: + raise ValueError('DiscreteMapping has no attached callable') + return self._callable_map + + def has_callable_mapping(self): + # WP07e: always True for any DiscreteMapping built through the public + # API (__new__ requires a callable, set_callable_mapping cannot detach + # it) -- exists so callers can ask without relying on + # get_callable_mapping()'s ValueError as control flow. + return self._callable_map is not None + + def set_callable_mapping(self, F): + # WP07d: unlike SymbolicMapping.set_callable_mapping (which this would + # otherwise inherit unguarded), a DiscreteMapping's callable is fixed at + # construction -- it is part of `_hashable_content`, so silently + # reattaching it would let `geo.mappings[name]` (read by assembly) and + # `self.get_callable_mapping()` (read by point-evaluation consumers) + # diverge without the object's identity (hash/equality) ever changing. + raise TypeError( + "DiscreteMapping's callable is fixed at construction (it is part " + "of its identity -- see _hashable_content) and cannot be " + "reattached; build a new DiscreteMapping(F, name=...) instead.") + + def copy(self): + # type(self)(name, ldim=...) -- SymbolicMapping.copy()'s call -- does not + # match DiscreteMapping.__new__, so rebuild explicitly (as + # InterfaceMapping / InverseMapping do), then carry the attributes + # SymbolicMapping.copy() preserves. + obj = DiscreteMapping(self._callable_map, self.name, + ldim=self.ldim, pdim=self.pdim) + obj._coordinates = self._coordinates + obj._logical_coordinates = self._logical_coordinates + obj._is_plus = self._is_plus + obj._is_minus = self._is_minus + return obj + + def _hashable_content(self): + # The wrapped callable is part of the identity: two same-named + # DiscreteMappings over different geometries must not compare equal + # (else MappedDomain's @cacheit conflates them). __new__ attaches + # `_callable_map` before any Indexed(self, i) / Jacobian(self) is built, + # so this is well-defined from the first hash. + cm = self._callable_map + key = cm if getattr(type(cm), '__hash__', None) else id(cm) + return super()._hashable_content() + (key,) + + @property + def func(self): + # Base `Basic.func` is `type(self)`, so `expr.func(*expr.args)` would be + # `DiscreteMapping(Symbol(name), Tuple(pdim))` -> TypeError. A + # DiscreteMapping has no symbolic sub-structure (the callable is not in + # `.args`), so the faithful rebuild IS the identity. sympy does this + # `func(*args)` rebuild while walking an expression that contains + # `self[i]` (e.g. `sympy.simplify` inside `Jacobian(M).inv()` when a + # bilinear form is lowered on a DiscreteMapping-carried domain). + return lambda *args, **kwargs: self + +#============================================================================== +class InverseMapping(StructuralMapping, metaclass=_MappingABCMeta): + """ Symbolic inverse of a mapping: F^{-1}. Not point-evaluable. """ + def __new__(cls, mapping): - assert isinstance(mapping, Mapping) - name = mapping.name - ldim = mapping.ldim - pdim = mapping.pdim - coords = mapping.logical_coordinates - jacobian = mapping.jacobian.inv() - return Mapping.__new__(cls, name, ldim=ldim, pdim=pdim, coordinates=coords, jacobian=jacobian) + assert isinstance(mapping, SymbolicMapping) + # Build the IndexedBase directly (WP06d-4a: no longer routing through + # Mapping.__new__, which wrapped `coordinates` Symbols in Symbol(...) + # again -- the pre-existing `Symbol(Symbol(...))` TypeError). + obj = IndexedBase.__new__(cls, mapping.name, shape=mapping.pdim) + lcoords = mapping._logical_coordinates # raw Tuple + obj._name = mapping.name + obj._ldim = mapping.ldim + obj._pdim = mapping.pdim + obj._coordinates = lcoords + obj._logical_coordinates = lcoords + obj._jacobian = mapping.jacobian_symbol.inv() + obj._is_minus = None + obj._is_plus = None + obj._base_mapping = mapping + cls._init_symbolic_jacobian(obj) + return obj + + @property + def ldim(self): + return self._ldim + + @property + def pdim(self): + return self._pdim + + @property + def jacobian_symbol(self): + return self._jacobian + + @property + def is_analytical(self): + return self._base_mapping.is_analytical + + def copy(self): + return InverseMapping(self._base_mapping) #============================================================================== class JacobianSymbol(MatrixSymbolicExpr): _axis = None def __new__(cls, mapping, axis=None): - assert isinstance(mapping, Mapping) + assert isinstance(mapping, SymbolicMapping) # incl. structural mappings (WP06d-4a) if axis is not None: assert isinstance(axis, (int, Integer)) obj = MatrixSymbolicExpr.__new__(cls, mapping) @@ -458,7 +1228,7 @@ def __hash__(self): return hash(self._hashable_content()) def _eval_subs(self, old, new): - if isinstance(new, Mapping): + if isinstance(new, SymbolicMapping): if self.axis is not None: obj = JacobianSymbol(new, self.axis) else: @@ -477,7 +1247,7 @@ class JacobianInverseSymbol(MatrixSymbolicExpr): _axis = None is_Matrix = False def __new__(cls, mapping, axis=None): - assert isinstance(mapping, Mapping) + assert isinstance(mapping, SymbolicMapping) # incl. structural mappings (WP06d-4a) if axis is not None: assert isinstance(axis, int) obj = MatrixSymbolicExpr.__new__(cls, mapping) @@ -509,7 +1279,7 @@ def _sympystr(self, printer): return 'Jacobian({})**(-1)'.format(sstr(self.mapping.name)) #============================================================================== -class InterfaceMapping(Mapping): +class InterfaceMapping(StructuralMapping, metaclass=_MappingABCMeta): """ InterfaceMapping is used to represent a mapping in the interface. @@ -522,18 +1292,32 @@ class InterfaceMapping(Mapping): """ def __new__(cls, minus, plus): - assert isinstance(minus, Mapping) - assert isinstance(plus, Mapping) + # Interface legs are patch mappings: an AnalyticMapping subclass, or a + # plain SymbolicMapping (e.g. built by Domain.from_file). The body needs + # copy() / set_plus_minus(), both on SymbolicMapping since WP06d-4c. + assert isinstance(minus, SymbolicMapping) + assert isinstance(plus, SymbolicMapping) minus = minus.copy() plus = plus.copy() minus.set_plus_minus(minus=True) plus.set_plus_minus(plus=True) - name = '{}|{}'.format(str(minus.name), str(plus.name)) - obj = Mapping.__new__(cls, name, ldim=minus.ldim, pdim=minus.pdim) - obj._minus = minus - obj._plus = plus + name = '{}|{}'.format(str(minus.name), str(plus.name)) + ldim, pdim = minus.ldim, minus.pdim + + # WP06d-4a: build the IndexedBase directly instead of via Mapping.__new__. + obj = IndexedBase.__new__(cls, name, shape=pdim) + obj._name = name + obj._ldim = ldim + obj._pdim = pdim + obj._coordinates = tuple(Symbol(u, real=True) for u in ['x', 'y', 'z'][:pdim]) + obj._logical_coordinates = Tuple(*(Symbol(u, real=True) for u in ['x1', 'x2', 'x3'][:ldim])) + obj._is_minus = None + obj._is_plus = None + obj._minus = minus + obj._plus = plus + cls._init_symbolic_jacobian(obj) return obj @property @@ -544,10 +1328,21 @@ def minus(self): def plus(self): return self._plus + @property + def ldim(self): + return self._ldim + + @property + def pdim(self): + return self._pdim + @property def is_analytical(self): return self.minus.is_analytical and self.plus.is_analytical + def copy(self): + return InterfaceMapping(self._minus, self._plus) + def _eval_subs(self, old, new): minus = self.minus.subs(old, new) plus = self.plus.subs(old, new) @@ -557,19 +1352,32 @@ def _eval_simplify(self, **kwargs): return self #============================================================================== -class MultiPatchMapping(Mapping): +class MultiPatchMapping(StructuralMapping, metaclass=_MappingABCMeta): def __new__(cls, dic): - assert isinstance( dic, dict) - return Basic.__new__(cls, dic) + if isinstance(dic, Dict): + dic = dict(dic.items()) + if not isinstance(dic, dict): + raise TypeError('Expecting a dictionary of patch mappings') + obj = _new_basic(cls, dic) + obj._mappings = dict(dic) + # WP06d-4a bug fix: MultiPatchMapping.__new__ never set _name, so any + # domain call or `==` that reached `name` raised AttributeError. + # An empty dict is degenerate but constructed fine pre-WP06d-4a. + first = next(iter(dic.values()), None) + obj._name = '|'.join(str(m.name) for m in dic.values()) + obj._coordinates = first._coordinates if first is not None else None + obj._logical_coordinates = first._logical_coordinates if first is not None else None + obj._is_minus = None + obj._is_plus = None + # NB: no _init_symbolic_jacobian -- MultiPatchMapping uses Basic.__new__ + # (not IndexedBase), so it can't be indexed; pre-WP06d-4a it never went + # through Mapping.__new__ either, so jacobian_expr stays None. + return obj @property def mappings(self): - return self.args[0] - - @property - def is_analytical(self): - return all(a.is_analytical for a in self.mappings.values()) + return self._mappings @property def ldim(self): @@ -583,14 +1391,17 @@ def pdim(self): def is_analytical(self): return all(e.is_analytical for e in self.mappings.values()) + def copy(self): + return MultiPatchMapping(dict(self.mappings)) + def _eval_subs(self, old, new): return self def _eval_simplify(self, **kwargs): return self - def __hash__(self): - return hash((*self.mappings.values(), *self.mappings.keys())) + def _hashable_content(self): + return self.args def _sympystr(self, printer): sstr = printer.doprint @@ -603,7 +1414,10 @@ class MappedDomain(BasicDomain): @cacheit def __new__(cls, mapping, logical_domain): - assert(isinstance(mapping, Mapping)) + # SymbolicMapping, not Mapping: since WP06d-4a the structural mappings + # (InterfaceMapping / MultiPatchMapping / InverseMapping) are callable + # on a domain but no longer subclass Mapping. + assert(isinstance(mapping, SymbolicMapping)) assert(isinstance(logical_domain, BasicDomain)) if isinstance(logical_domain, Domain): kwargs = dict( @@ -701,7 +1515,7 @@ def __new__(cls, u, mapping=None): else: mapping = space.domain.mapping - J = mapping.jacobian + J = mapping.jacobian_symbol if isinstance(kind, (UndefinedSpaceType, H1SpaceType)): expr = el @@ -762,8 +1576,12 @@ def eval(cls, F): the jacobian matrix """ - if not isinstance(F, Mapping): - raise TypeError('> Expecting a Mapping object') + # SymbolicMapping, not Mapping: WP06d-4a's structural mappings ask for + # their own Jacobian(self) at construction (via + # StructuralMapping._init_symbolic_jacobian), and they are no longer + # Mapping subclasses. + if not isinstance(F, SymbolicMapping): + raise TypeError('> Expecting a SymbolicMapping object') if F.jacobian_expr is not None: return F.jacobian_expr @@ -862,8 +1680,11 @@ def eval(cls, F, v): the contravariant transformation """ - if not isinstance(F, Mapping): - raise TypeError('> Expecting a Mapping') + # SymbolicMapping, not Mapping: consistent with Covariant.eval (no + # guard) and Jacobian.eval -- WP06d-4a's structural mappings are valid + # here (Hdiv push-forward over a multipatch/broken domain). + if not isinstance(F, SymbolicMapping): + raise TypeError('> Expecting a SymbolicMapping') if not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)): raise TypeError('> Expecting a tuple, list, Tuple, Matrix') @@ -1009,7 +1830,7 @@ def eval(cls, expr, domain, **options): else: arg = cls.eval(arg, domain) - return mapping.jacobian.inv().T*grad(arg) + return mapping.jacobian_symbol.inv().T*grad(arg) elif isinstance(expr, curl): arg = expr.args[0] @@ -1029,7 +1850,7 @@ def eval(cls, expr, domain, **options): arg = cls.eval(arg, domain) if isinstance(arg, PullBack) and isinstance(arg.kind, HcurlSpaceType): - J = mapping.jacobian + J = mapping.jacobian_symbol arg = arg.test if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)): arg = type(expr.args[0])(arg) @@ -1059,20 +1880,20 @@ def eval(cls, expr, domain, **options): arg = cls.eval(arg, domain) if isinstance(arg, PullBack) and isinstance(arg.kind, HdivSpaceType): - J = mapping.jacobian + J = mapping.jacobian_symbol arg = arg.test if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)): arg = type(expr.args[0])(arg) return (1/J.det())*div(arg) elif isinstance(arg, PullBack): - return SymbolicTrace(mapping.jacobian.inv().T*grad(arg.test)) + return SymbolicTrace(mapping.jacobian_symbol.inv().T*grad(arg.test)) else: raise NotImplementedError('TODO') elif isinstance(expr, laplace): arg = expr.args[0] v = cls.eval(grad(arg), domain) - v = mapping.jacobian.inv().T*grad(v) + v = mapping.jacobian_symbol.inv().T*grad(v) return SymbolicTrace(v) # elif isinstance(expr, hessian): @@ -1220,7 +2041,7 @@ def eval(cls, expr, domain, **options): assert domain is not None if expr.is_domain_integral: - J = mapping.jacobian + J = mapping.jacobian_symbol det = sqrt((J.T*J).det()) else: axis = domain.axis @@ -1253,7 +2074,7 @@ def eval(cls, expr, domain, **options): elif isinstance(expr, DomainExpression): domain = expr.target - J = domain.mapping.jacobian + J = domain.mapping.jacobian_symbol newexpr = cls.eval(expr.expr, domain) newexpr = TerminalExpr(newexpr, domain=domain) domain = domain.logical_domain @@ -1355,7 +2176,9 @@ def eval(cls, *_args, **kwargs): elif isinstance(expr, Indexed): base = expr.base - if isinstance(base, Mapping): + # SymbolicMapping, not Mapping: WP06d-4a's structural mappings are + # IndexedBase-backed but no longer Mapping (mirrors the L1800 check). + if isinstance(base, SymbolicMapping): if expr.indices[0] == 0: name = 'x' elif expr.indices[0] == 1: @@ -1400,7 +2223,7 @@ def eval(cls, *_args, **kwargs): code += k*n return cls.eval(atom, code=code) - elif isinstance(expr, Mapping): + elif isinstance(expr, SymbolicMapping): return Symbol(expr.name) # ... this must be done here, otherwise codegen for FEM will not work diff --git a/sympde/topology/space.py b/sympde/topology/space.py index cce1f556..6495bb6a 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -16,6 +16,7 @@ from sympde.core.utils import random_string from sympde.core.basic import CalculusFunction from sympde.core.basic import _coeffs_registery +from sympde.core.basic import _new_basic from .basic import BasicDomain, Union, Interval from .datatype import SpaceType, dtype_space_registry @@ -219,6 +220,10 @@ def __mul__(self, other): def __hash__(self): return hash((self.name, self.domain, self.shape, self.kind)) + + def _hashable_content(self): + return (self.name, self.domain, self.shape, self.kind) + #============================================================================== class ScalarFunctionSpace(BasicFunctionSpace): """ @@ -438,19 +443,14 @@ class ScalarFunction(Symbol): def __new__(cls, space, name): if not isinstance(space, ScalarFunctionSpace): raise ValueError('Expecting a ScalarFunctionSpace') - obj = Expr.__new__(cls) + obj = Symbol.__xnew__(cls, name, commutative=True) obj._space = space - obj._name = name return obj @property def space(self): return self._space - @property - def name(self): - return self._name - @property def ldim(self): return self.space.ldim @@ -469,8 +469,9 @@ def _sympystr(self, printer): sstr = printer.doprint return sstr(self.name) - def __hash__(self): - return hash((self.name, self.space)) + def _hashable_content(self): + return Symbol._hashable_content(self) + (self.space,) + #============================================================================== # this class is needed, otherwise sympy will convert VectorFunction to # IndexedBase @@ -547,19 +548,14 @@ class VectorFunction(Symbol, IndexedBase): def __new__(cls, space, name): if not isinstance(space, VectorFunctionSpace): raise ValueError('Expecting a VectorFunctionSpace') - obj = Expr.__new__(cls) + obj = Symbol.__xnew__(cls, name, commutative=False) obj._space = space - obj._name = name return obj @property def space(self): return self._space - @property - def name(self): - return self._name - @property def shape(self): # we return a list to make it compatible with IndexedBase sympy object @@ -595,8 +591,8 @@ def _sympystr(self, printer): sstr = printer.doprint return sstr(self.name) - def __hash__(self): - return hash((self.name, self.space)) + def _hashable_content(self): + return Symbol._hashable_content(self) + (self.space,) #============================================================================== # this is implemented as a function, it would be better to have it as a class @@ -646,7 +642,7 @@ def __new__(cls, expr, boundary, order=0, **options): if options.pop('evaluate',True): return cls.eval(expr, boundary, order) - obj = Basic.__new__(cls, expr, boundary, order) + obj = _new_basic(cls, expr, boundary, order) obj.is_commutative = expr.is_commutative return obj @@ -660,7 +656,7 @@ def boundary(self): @property def order(self): - return self._args[2] + return int(self._args[2]) @classmethod def eval(cls, expr, boundary, order): diff --git a/sympde/topology/tests/test_callable_mapping.py b/sympde/topology/tests/test_callable_mapping.py index 23462c3d..41a20fb1 100644 --- a/sympde/topology/tests/test_callable_mapping.py +++ b/sympde/topology/tests/test_callable_mapping.py @@ -1,6 +1,7 @@ import numpy as np +import pytest -from sympde.topology.mapping import Mapping, BasicCallableMapping +from sympde.topology.mapping import SymbolicMapping, BasicCallableMapping from sympde.topology.analytical_mapping import IdentityMapping, AffineMapping from sympde.topology.analytical_mapping import PolarMapping @@ -506,9 +507,12 @@ def ldim(self): def pdim(self): return self._ndim - F = Mapping('F', ldim = 3, pdim = 3) # Declare undefined symbolic mapping + F = SymbolicMapping('F', ldim = 3, pdim = 3) # Declare undefined symbolic mapping f = UserIdentity(3) # Create user-defined callable mapping - F.set_callable_mapping(f) # Attach callable mapping to symbolic mapping + # set_callable_mapping is deprecated; DiscreteMapping(f, 'F') is the + # current route to a point-evaluable symbolic identity. + with pytest.warns(DeprecationWarning): + F.set_callable_mapping(f) # Attach callable mapping to symbolic mapping assert F.get_callable_mapping() is f assert f(4, 5, 6) == (4, 5, 6) diff --git a/sympde/topology/tests/test_derivatives.py b/sympde/topology/tests/test_derivatives.py index 3694adfe..e3963fdc 100644 --- a/sympde/topology/tests/test_derivatives.py +++ b/sympde/topology/tests/test_derivatives.py @@ -12,7 +12,7 @@ from sympde.topology import get_max_partial_derivatives from sympde.topology import ScalarFunctionSpace from sympde.topology import (dx, dy, dz) -from sympde.topology import Mapping +from sympde.topology import SymbolicMapping def indices_as_str(a): @@ -30,7 +30,7 @@ def test_partial_derivatives_1(): # ... domain = Domain('Omega', dim=2) - M = Mapping('M', dim=2) + M = SymbolicMapping('M', dim=2) mapped_domain = M(domain) @@ -78,7 +78,7 @@ def test_partial_derivatives_2(): # ... domain = Domain('Omega', dim=2) - M = Mapping('M', dim=2) + M = SymbolicMapping('M', dim=2) mapped_domain = M(domain) @@ -130,6 +130,99 @@ def test_partial_derivatives_2(): # ... +# ... +def test_logical_derivative_through_symbolic_mapping_index(): + # /code-review finding 1 (WP06d-4c-1a): `_DifferentialOperator.eval`'s + # logical chain-rule branch matched `expr.atoms(Mapping)`, which no longer + # catches a bare `SymbolicMapping` (the constructor WP06d-4c recommends and + # migrated every test to). The chain-rule term through the mapping component + # `M[i]` was then silently dropped and `dx1(M[0]**2)` collapsed to 0. + from sympde.topology import dx1, dx2 + + M = SymbolicMapping('M', dim=2) + + assert dx1(M[0]**2) == 2 * M[0] * dx1(M[0]) + assert dx2(M[1]**2) == 2 * M[1] * dx2(M[1]) + + N = SymbolicMapping('N', dim=3) + assert dx1(N[2]**3) == 3 * N[2]**2 * dx1(N[2]) +# ... + + +# ... +def test_logical_derivative_through_interface_mapping_component(): + # /code-review finding 1 (WP06d-4c-1b): the 06d-4c-1a widening to + # `expr.atoms(SymbolicMapping)` now also selects structural mappings. The + # indexable ones (InterfaceMapping, InverseMapping) must still be + # chain-ruled through -- the `_shape` filter added in 4c-1b keeps them. + from sympde.topology import dx1, InterfaceMapping, IdentityMapping + + itf = InterfaceMapping(IdentityMapping('A', dim=2), IdentityMapping('B', dim=2)) + assert dx1(itf[0]**2) == 2 * itf[0] * dx1(itf[0]) + + +def test_multipatch_mapping_is_excluded_from_the_chain_rule_branch(): + # WP06d-4c-1b: the chain-rule branch indexes the mapping (`M[i]`), so it + # filters `expr.atoms(SymbolicMapping)` to indexable mappings via + # `getattr(m, '_shape', None) is not None`. MultiPatchMapping is built via + # Basic.__new__ and has no `_shape`; without the filter `M[i]` raised + # AttributeError. (It cannot actually be built into a scalar Expr reaching + # that branch -- its `.args` is a raw dict, which breaks Expr.is_number + # first -- so the filter is a defensive guard; lock the discriminator here.) + from sympde.topology import (MultiPatchMapping, InterfaceMapping, + IdentityMapping) + + mp = MultiPatchMapping({'p': IdentityMapping('F', dim=2)}) + itf = InterfaceMapping(IdentityMapping('A', dim=2), IdentityMapping('B', dim=2)) + assert getattr(mp, '_shape', None) is None + assert getattr(itf, '_shape', None) is not None +# ... + + +# ... +def test_logical_derivative_sums_over_all_mappings(): + # D-1: `_DifferentialOperator.eval`'s logical chain-rule branch used to pick + # ONE mapping out of `expr.atoms(SymbolicMapping)` and treat the rest as + # constants -- a silently wrong derivative when two or more mappings' indexed + # components appear in the same expression. + from sympy import sqrt, simplify + from sympde.topology import dx1 + + F = SymbolicMapping('F', dim=2) + G = SymbolicMapping('G', dim=2) + + r = dx1(sqrt(F[0]**2 + G[0]**2)) + expected = (F[0]*dx1(F[0]) + G[0]*dx1(G[0])) / sqrt(F[0]**2 + G[0]**2) + assert simplify(r - expected) == 0 + assert dx1(F[0]) in r.atoms(type(dx1(F[0]))) + assert dx1(G[0]) in r.atoms(type(dx1(G[0]))) + + +def test_logical_derivative_single_mapping_is_unchanged(): + # D-1 generalisation must not perturb the common one-mapping path. + from sympy import sqrt, simplify + from sympde.topology import dx1, dx2 + + F = SymbolicMapping('F', dim=2) + + r = dx1(sqrt(F[0]**2 + F[1]**2)) + assert simplify(r - (F[0]*dx1(F[0]) + F[1]*dx1(F[1])) / sqrt(F[0]**2 + F[1]**2)) == 0 + assert dx2(F[0]**3) == 3 * F[0]**2 * dx2(F[0]) + + +def test_logical_derivative_through_patch_map_and_interface_mapping(): + # D-1: a patch SymbolicMapping and an InterfaceMapping routinely co-occur; + # both must contribute a chain-rule term. + from sympde.topology import dx1, InterfaceMapping, IdentityMapping + + F = SymbolicMapping('F', dim=2) + itf = InterfaceMapping(IdentityMapping('A', dim=2), IdentityMapping('B', dim=2)) + + r = dx1(F[0]*itf[0]) + assert r == F[0]*dx1(itf[0]) + itf[0]*dx1(F[0]) +# ... + + #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== diff --git a/sympde/topology/tests/test_logical_expr.py b/sympde/topology/tests/test_logical_expr.py index b66d3d86..e7801ad5 100644 --- a/sympde/topology/tests/test_logical_expr.py +++ b/sympde/topology/tests/test_logical_expr.py @@ -9,7 +9,7 @@ from sympde.calculus import grad, dot, inner, rot, div from sympde.calculus import laplace, bracket, convect from sympde.calculus import jump, avg, Dn, minus, plus -from sympde.topology import Domain, Mapping, Square +from sympde.topology import Domain, SymbolicMapping, Square from sympde.topology import dx, dy from sympde.topology import dx1, dx2, dx3 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace @@ -23,6 +23,7 @@ from sympde.topology import CollelaMapping2D from sympde.topology import TorusMapping from sympde.topology import TwistedTargetMapping +from sympde.topology.mapping import get_logical_test_function from sympde.expr import BilinearForm, integral from sympde.calculus import grad, div, curl, dot @@ -33,7 +34,7 @@ def test_logical_expr_1d_1(): dim = 1 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) domain = M(Domain('Omega', dim=dim)) alpha = Constant('alpha') @@ -80,7 +81,7 @@ def test_logical_expr_1d_1(): def test_symbolic_expr_1d_1(): dim = 1 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) domain = M(Domain('Omega', dim=dim)) V = ScalarFunctionSpace('V', domain, kind='h1') @@ -149,7 +150,7 @@ def test_symbolic_expr_1d_1(): def test_logical_expr_2d_1(): dim = 2 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) domain = M(Domain('Omega', dim=dim)) alpha = Constant('alpha') @@ -215,8 +216,8 @@ def test_logical_expr_2d_2(): A = Square('A') B = Square('B') - M1 = Mapping('M1', dim=dim) - M2 = Mapping('M2', dim=dim) + M1 = SymbolicMapping('M1', dim=dim) + M2 = SymbolicMapping('M2', dim=dim) D1 = M1(A) D2 = M2(B) @@ -268,8 +269,8 @@ def test_logical_expr_2d_3(): A = Square('A') B = Square('B') - M1 = Mapping('M1', dim=dim) - M2 = Mapping('M2', dim=dim) + M1 = SymbolicMapping('M1', dim=dim) + M2 = SymbolicMapping('M2', dim=dim) D1 = M1(A) D2 = M2(B) @@ -293,7 +294,7 @@ def test_logical_expr_2d_3(): #============================================================================== def test_symbolic_expr_2d_1(): dim = 2 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) domain = M(Domain('Omega', dim=dim)) V = ScalarFunctionSpace('V', domain, kind='h1') @@ -369,7 +370,7 @@ def test_symbolic_expr_2d_1(): def test_logical_expr_3d_1(): dim = 3 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) domain = M(Domain('Omega', dim=dim)) alpha = Constant('alpha') @@ -422,17 +423,16 @@ def test_logical_expr_3d_2(): dim = 3 domain = Domain('Omega', dim=dim) - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) mapped_domain = M(domain) - V = ScalarFunctionSpace('V' , domain, kind='h1') VM = ScalarFunctionSpace('VM', mapped_domain, kind='h1') - u,v = elements_of(V, names='u,v') um,vm = elements_of(VM, names='u,v') + u, v = [get_logical_test_function(field) for field in (um, vm)] - J = M.jacobian + J = M.jacobian_symbol a = dot(grad(um),grad(vm)) e = LogicalExpr(a, mapped_domain) @@ -444,17 +444,16 @@ def test_logical_expr_3d_3(): dim = 3 domain = Domain('Omega', dim=dim) - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) mapped_domain = M(domain) - V = VectorFunctionSpace('V' , domain, kind='hcurl') VM = VectorFunctionSpace('VM', mapped_domain, kind='hcurl') - u,v = elements_of(V, names='u,v') um,vm = elements_of(VM, names='u,v') + u, v = [get_logical_test_function(field) for field in (um, vm)] - J = M.jacobian + J = M.jacobian_symbol a = dot(curl(um), curl(vm)) e = LogicalExpr(a, mapped_domain) @@ -465,17 +464,16 @@ def test_logical_expr_3d_4(): dim = 3 domain = Domain('Omega', dim=dim) - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) mapped_domain = M(domain) - V = VectorFunctionSpace('V' , domain, kind='hdiv') VM = VectorFunctionSpace('VM', mapped_domain, kind='hdiv') - u,v = elements_of(V, names='u,v') um,vm = elements_of(VM, names='u,v') + u, v = [get_logical_test_function(field) for field in (um, vm)] - J = M.jacobian + J = M.jacobian_symbol a = div(um)*div(vm) e = LogicalExpr(a, mapped_domain) @@ -486,16 +484,15 @@ def test_logical_expr_3d_5(): dim = 3 domain = Domain('Omega', dim=dim) - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) mapped_domain = M(domain) - V = VectorFunctionSpace('V' , domain, kind='hcurl') VM = VectorFunctionSpace('VM', mapped_domain, kind='hcurl') - J = M.jacobian - u,v = elements_of(V, names='u,v') + J = M.jacobian_symbol um,vm = elements_of(VM, names='u,v') + u, v = [get_logical_test_function(field) for field in (um, vm)] int_md = lambda expr: integral(mapped_domain , expr) int_ld = lambda expr: integral(domain , expr) @@ -510,7 +507,7 @@ def test_logical_expr_3d_5(): #============================================================================== def test_symbolic_expr_3d_1(): dim = 3 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) domain = M(Domain('Omega', dim=dim)) V = ScalarFunctionSpace('V', domain, kind='h1') @@ -617,10 +614,12 @@ def test_identity_mapping_2d_2(): V = ScalarFunctionSpace('V', domain, kind='h1') u = element_of(V, name='u') + logical_u = get_logical_test_function(u) # ... - assert(LogicalExpr(dx(u), domain) == dx1(u)) - assert(LogicalExpr(dy(u), domain) == dx2(u)) + assert logical_u.space.domain == domain.logical_domain + assert(LogicalExpr(dx(u), domain) == dx1(logical_u)) + assert(LogicalExpr(dy(u), domain) == dx2(logical_u)) # ... #============================================================================== diff --git a/sympde/topology/tests/test_mapping.py b/sympde/topology/tests/test_mapping.py index 8d7c6798..fbd29fcb 100644 --- a/sympde/topology/tests/test_mapping.py +++ b/sympde/topology/tests/test_mapping.py @@ -7,19 +7,38 @@ from sympy.tensor import IndexedBase from sympy import symbols, simplify -from sympde.topology import Mapping, MappedDomain, AffineMapping +from sympde.topology import IdentityMapping, SymbolicMapping, MappedDomain, AffineMapping from sympde.topology import dx, dy, dz from sympde.topology import dx1, dx2, dx3 -from sympde.topology import Domain +from sympde.topology import Domain, Square -from sympde.topology.mapping import Jacobian, Covariant, Contravariant +from sympde.topology.mapping import ( + Contravariant, + Covariant, + Jacobian, + MultiPatchMapping, +) + + +def test_multipatch_mapping_reconstruction(): + patch = Square('A') + mapping = IdentityMapping('F', dim=2) + replacement = IdentityMapping('G', dim=2) + multipatch_mapping = MultiPatchMapping({patch: mapping}) + + rebuilt = multipatch_mapping.func(*multipatch_mapping.args) + replaced = multipatch_mapping.xreplace({mapping: replacement}) + + assert rebuilt == multipatch_mapping + assert rebuilt.mappings == multipatch_mapping.mappings + assert replaced.mappings == {patch: replacement} # ... def test_mapping_1d(): print('============ test_mapping_1d ==============') dim = 1 - F = Mapping('F', dim=dim) + F = SymbolicMapping('F', dim=dim) assert(F.name == 'F') @@ -40,7 +59,7 @@ def test_mapping_2d(): dim = 2 - F = Mapping('F', dim=dim) + F = SymbolicMapping('F', dim=dim) a,b = symbols('a b') ab = Tuple(a, b) @@ -80,7 +99,7 @@ def test_mapping_3d(): dim = 3 - F = Mapping('F', dim=dim) + F = SymbolicMapping('F', dim=dim) a,b,c = symbols('a b c') abc = Tuple(a, b, c) @@ -125,7 +144,7 @@ def test_mapping_2d_2(): print('============ test_mapping_2d_2 ==============') dim = 2 - F = Mapping('F', dim=dim) + F = SymbolicMapping('F', dim=dim) domain = Domain('Omega', dim=dim) D = F(domain) @@ -157,4 +176,3 @@ def teardown_module(): def teardown_function(): from sympy.core import cache cache.clear_cache() - diff --git a/sympde/topology/tests/test_mapping_hierarchy.py b/sympde/topology/tests/test_mapping_hierarchy.py new file mode 100644 index 00000000..b7336896 --- /dev/null +++ b/sympde/topology/tests/test_mapping_hierarchy.py @@ -0,0 +1,1064 @@ +# coding: utf-8 +""" +Work-package 01: the new mapping base classes (SymbolicMapping, DefinedMapping, +StructuralMapping) introduced alongside the legacy hierarchy. + +These checks lock the DefinedMapping interface and the abstractness of +DefinedMapping / StructuralMapping without exercising sympy object construction +(that is covered by later work-packages). +""" +import inspect +import warnings + +import pytest + +from sympde.topology import SymbolicMapping, DefinedMapping, StructuralMapping +from sympde.topology import mapping as _mapping_mod +from sympde.topology.mapping import BasicCallableMapping + + +def test_new_classes_are_exported(): + # the star-import from the package must expose the very objects defined in + # the module, and the module must advertise them in __all__ + assert SymbolicMapping is _mapping_mod.SymbolicMapping + assert DefinedMapping is _mapping_mod.DefinedMapping + assert StructuralMapping is _mapping_mod.StructuralMapping + for name in ('SymbolicMapping', 'DefinedMapping', 'StructuralMapping'): + assert name in _mapping_mod.__all__ + + +def test_hierarchy_links(): + assert issubclass(DefinedMapping, SymbolicMapping) + assert issubclass(StructuralMapping, SymbolicMapping) + + +def test_defined_mapping_interface_covers_basic_callable_mapping(): + # DefinedMapping supersedes BasicCallableMapping: every name of the old + # interface is still exposed. 06d-4c: ldim / pdim / __call__ are now + # provided concretely by SymbolicMapping, so only the numeric point-eval + # methods stay abstract until a concrete subclass implements them. + for name in BasicCallableMapping.__abstractmethods__: + assert hasattr(DefinedMapping, name) + assert {'jacobian', 'jacobian_inv', 'metric', 'metric_det'}.issubset( + DefinedMapping.__abstractmethods__) + + +def test_symbolic_mapping_is_concrete(): + assert not inspect.isabstract(SymbolicMapping) + + +@pytest.mark.parametrize('cls', [DefinedMapping, StructuralMapping]) +def test_abstract_classes_reject_instantiation(cls): + assert inspect.isabstract(cls) + with pytest.raises(TypeError, match='abstract'): + cls('F') + + +def test_concrete_subclasses_lose_abstractness_without_instantiation(): + class ConcreteDefined(DefinedMapping): + def __call__(self, *eta): return eta + def jacobian(self, *eta): return None + def jacobian_inv(self, *eta): return None + def metric(self, *eta): return None + def metric_det(self, *eta): return None + ldim = 2 + pdim = 2 + + class ConcreteStructural(StructuralMapping): + ldim = 2 + pdim = 2 + + assert not inspect.isabstract(ConcreteDefined) + assert not inspect.isabstract(ConcreteStructural) + assert ConcreteDefined.__abstractmethods__ == frozenset() + assert ConcreteStructural.__abstractmethods__ == frozenset() + + +# -- work-package 02a: the symbolic `jacobian` property moved to `jacobian_symbol` + +def test_jacobian_symbol_holds_the_symbolic_jacobian(): + from sympde.topology import IdentityMapping + F = IdentityMapping('F', dim=2) + assert F.jacobian_symbol is F._jacobian + assert type(F.jacobian_symbol).__name__ == 'JacobianSymbol' + + +def test_legacy_jacobian_property_is_a_deprecated_alias(): + from sympde.topology import Mapping + # `.jacobian` (the deprecated symbolic property) lives on `Mapping` only -- + # `SymbolicMapping` deliberately does not carry it. Since WP D8, + # `AnalyticMapping` no longer inherits `Mapping`, so nothing is "shadowed" + # any more: `.jacobian` is simply confined to `Mapping` and its legacy + # subclasses, while `AnalyticMapping.jacobian` is a distinct, unrelated + # numeric method (see test_analytic_mapping_is_directly_point_evaluable). + with pytest.warns(DeprecationWarning, match='SymbolicMapping'): # 06d-4c: ctor + F = Mapping('F', dim=2) + with pytest.warns(DeprecationWarning, match='jacobian_symbol'): + legacy = F.jacobian + assert legacy is F.jacobian_symbol + + +# -- work-package 02b: AnalyticMapping is a concrete, point-evaluable DefinedMapping + +def test_defined_mapping_inherits_basic_callable_mapping(): + # D5: the point-eval interface is now stated once, on BasicCallableMapping. + # 06d-4c: SymbolicMapping provides ldim / pdim / __call__ concretely, so + # DefinedMapping only still-abstracts the numeric point-eval methods. + assert issubclass(DefinedMapping, BasicCallableMapping) + assert DefinedMapping.__abstractmethods__ == frozenset( + {'jacobian', 'jacobian_inv', 'metric', 'metric_det'}) + + +def test_analytic_mapping_hierarchy(): + from sympde.topology import AnalyticMapping, Mapping + # WP D8: AnalyticMapping no longer inherits the deprecated Mapping. + assert not issubclass(AnalyticMapping, Mapping) + assert issubclass(AnalyticMapping, DefinedMapping) + assert issubclass(AnalyticMapping, BasicCallableMapping) + assert not inspect.isabstract(AnalyticMapping) + + +def test_analytical_gallery_reparented_onto_analytic_mapping(): + from sympde.topology import AnalyticMapping, Mapping + from sympde.topology import (IdentityMapping, AffineMapping, PolarMapping, + TargetMapping, CzarnyMapping, CollelaMapping2D, + TorusMapping) + for cls in (IdentityMapping, AffineMapping, PolarMapping, TargetMapping, + CzarnyMapping, CollelaMapping2D, TorusMapping): + assert issubclass(cls, AnalyticMapping) + # WP D8: no longer under the deprecated Mapping; use SymbolicMapping + # as the "any symbolic mapping" check (same as 06d-4a for structural). + assert not issubclass(cls, Mapping) + assert issubclass(cls, SymbolicMapping) + + +# -- D8: AnalyticMapping no longer inherits the deprecated Mapping + +def test_analytic_mapping_no_longer_inherits_deprecated_mapping(): + from sympde.topology import AnalyticMapping, Mapping + from sympde.topology import PolarMapping + assert not issubclass(AnalyticMapping, Mapping) + assert Mapping not in PolarMapping.__mro__ + assert issubclass(AnalyticMapping, DefinedMapping) + assert issubclass(AnalyticMapping, SymbolicMapping) + assert not inspect.isabstract(AnalyticMapping) + # no C3 duplication of SymbolicMapping in the MRO + assert [c.__name__ for c in AnalyticMapping.__mro__].count('SymbolicMapping') == 1 + assert issubclass(Mapping, SymbolicMapping) + + +def test_analytic_expression_machinery_survives_the_move(): + from sympy import ImmutableDenseMatrix, eye + from sympde.topology import IdentityMapping, PolarMapping, TorusSurfaceMapping + + F = IdentityMapping('F', dim=2) + assert F.jacobian_expr == ImmutableDenseMatrix(eye(2)) + assert F.jacobian_inv_expr == ImmutableDenseMatrix(eye(2)) + assert F.metric_det_expr == 1 + + # no numeric constants given: all four stay symbolic Constants + from sympy import Symbol, cos, sin + from sympde.core.basic import Constant + P0 = PolarMapping('P0', dim=2) + assert set(a.name for a in P0.constants) == {'c1', 'c2', 'rmin', 'rmax'} + c1, c2, rmin, rmax = (Constant(n) for n in ('c1', 'c2', 'rmin', 'rmax')) + # the logical coordinates on the mapping are real-valued Symbols + x1, x2 = (Symbol(n, real=True) for n in ('x1', 'x2')) + expected = (c1 + (rmin*(1 - x1) + rmax*x1)*cos(x2), + c2 + (rmin*(1 - x1) + rmax*x1)*sin(x2)) + assert P0.expressions == expected + + # pdim != ldim: the _inv_jac branch is None + T = TorusSurfaceMapping('T', ldim=2, pdim=3, R0=1., a=0.3) + assert T.jacobian_expr.shape == (3, 2) + assert T.jacobian_inv_expr is None + assert T.metric_expr.shape == (2, 2) + + +def test_legacy_mapping_subclass_still_expands_expressions(): + from sympy import ImmutableDenseMatrix + from sympde.topology import Mapping, AnalyticMapping + + class M(Mapping): + _expressions = {'x': '2*x1', 'y': '3*x2'} + _ldim = 2 + _pdim = 2 + + with warnings.catch_warnings(): + warnings.simplefilter('error', DeprecationWarning) + m = M('m') # `cls is Mapping` gate: subclass construction is quiet + + assert m.jacobian_expr == ImmutableDenseMatrix([[2, 0], [0, 3]]) + assert m.metric_det_expr == 36 + assert not isinstance(m, AnalyticMapping) + + +def test_analytic_mapping_construction_is_warning_free(): + from sympde.topology import IdentityMapping, PolarMapping, TorusSurfaceMapping + with warnings.catch_warnings(): + warnings.simplefilter('error', DeprecationWarning) + F = IdentityMapping('F', dim=2) + PolarMapping('P', dim=2, c1=0., c2=0., rmin=.3, rmax=1.) + TorusSurfaceMapping('T', ldim=2, pdim=3, R0=1., a=0.3) + F.copy() + F.func(*F.args) + + +def test_analytic_mapping_is_directly_point_evaluable(): + from sympde.topology import IdentityMapping + F = IdentityMapping('F', dim=2) + cm = F.get_callable_mapping() + # point evaluation goes through the same callable mapping as before + assert F(0.3, 0.4) == cm(0.3, 0.4) + J = F.jacobian(0.3, 0.4) + assert type(J).__name__ != 'JacobianSymbol' # numeric, not symbolic + assert list(map(list, J)) == list(map(list, cm.jacobian(0.3, 0.4))) + assert F.metric_det(0.3, 0.4) == cm.metric_det(0.3, 0.4) + + +def test_analytic_mapping_still_callable_on_a_domain(): + from sympde.topology import IdentityMapping, Square + F = IdentityMapping('F', dim=2) + mapped = F(Square('D')) + assert type(mapped).__name__ in ('Domain', 'MappedDomain') + + +# -- D11: AnalyticMapping requires _expressions + +@pytest.mark.parametrize('kwargs', [ + dict(dim=2), + dict(ldim=2, pdim=3), + dict(ldim=2, pdim=2, evaluate=False), +]) +def test_bare_analytic_mapping_is_rejected_at_construction(kwargs): + # D11: a class with no `_expressions` cannot be point-evaluated, so + # constructing it is rejected outright, unconditionally on `evaluate` + # (the third case) -- the degenerate object must never come to exist. + from sympde.topology import AnalyticMapping + with pytest.raises(TypeError, match='analytical expressions') as exc_info: + AnalyticMapping('F', **kwargs) + assert 'SymbolicMapping' in str(exc_info.value) + + +def test_analytic_intermediate_base_without_expressions_is_definable_not_instantiable(): + # The guard lives in __new__, not __init_subclass__: an intermediate base + # without `_expressions` stays a legal *definition*, only instantiating it + # raises. + from sympde.topology import AnalyticMapping + + class Base(AnalyticMapping): + def helper(self): + return 'ok' + + with pytest.raises(TypeError, match='analytical expressions'): + Base('b', dim=2) + + class Leaf(Base): + _expressions = {'x': '2*x1', 'y': '3*x2'} + _ldim = 2 + _pdim = 2 + + L = Leaf('l') + assert L(1.0, 1.0) == (2.0, 3.0) + + +def test_construction_guard_spares_rebuilds_of_analytic_mappings(): + # R4: rebuild paths only rebuild instances of classes that already passed + # the guard, so they must not raise. Rebuild EQUALITY is deliberately not + # asserted here: P.func(*P.args) == P and a pickle round trip == P are + # already False at baseline for gallery mappings (pre-existing, numeric + # constants come back symbolic) -- out of scope for D11. + import pickle # nosec B403 -- only round-trips objects built by the test + from sympde.topology import PolarMapping, InterfaceMapping, IdentityMapping + + with warnings.catch_warnings(): + warnings.simplefilter('error', DeprecationWarning) + P = PolarMapping('P', dim=2, c1=0., c2=0., rmin=.3, rmax=1.) + + for rebuilt in (P.copy(), P.func(*P.args), pickle.loads(pickle.dumps(P))): # nosec B301 + assert isinstance(rebuilt, PolarMapping) + assert rebuilt.is_analytical is True + + itf = InterfaceMapping(IdentityMapping('A', dim=2), IdentityMapping('B', dim=2)) + assert isinstance(itf.minus, IdentityMapping) + + +def test_gallery_analytic_mappings_all_carry_expressions(): + # Filter by __module__, not AnalyticMapping.__subclasses__(): the latter + # would also pick up test-local subclasses (e.g. Base above), polluting + # the list depending on test order and xdist scheduling. + import sympde.topology.analytical_mapping as gallery + from sympde.topology import AnalyticMapping + + classes = [c for c in vars(gallery).values() + if isinstance(c, type) and issubclass(c, AnalyticMapping) + and c.__module__ == 'sympde.topology.analytical_mapping'] + assert len(classes) == 11 + for c in classes: + assert c._expressions is not None + + +# -- work-package 02b: post-review fixes + +def test_analytic_mapping_jacobian_still_shadows_symbolic_alias_but_jacobian_symbol_is_unaffected(): + # code-review fix #1, corrected: AnalyticMapping.jacobian (the numeric + # method) shadows the deprecated symbolic Mapping.jacobian property with no + # warning. But investigating the flagged psydac call sites + # (compute_boundary_jacobian / compute_normal_vector) showed + # SymbolicExpr() has never supported translating a JacobianSymbol -- even + # on a bare, pre-refactor Mapping, SymbolicExpr(mapping.jacobian) already + # raised NotImplementedError. So this is not a regression: those two + # psydac functions were already non-functional dead code (zero callers in + # the workspace) before this refactor. What we *can* guarantee here is + # that jacobian_symbol keeps returning the same JacobianSymbol object + # regardless of shadowing, so the eventual WP05 rename is behaviour + # preserving. + from sympde.topology import PolarMapping + F = PolarMapping('F', rmin=0, rmax=1, c1=0, c2=0) + assert type(F.jacobian_symbol).__name__ == 'JacobianSymbol' + assert type(F.jacobian) is not type(F.jacobian_symbol) # method vs symbol + + +def test_interface_mapping_copy_preserves_analytic_mapping_type(): + # code-review fix #2: Mapping.copy() used to hardcode Mapping(...), so + # InterfaceMapping (which copies both legs) silently downgraded an + # AnalyticMapping subclass to a plain, non-point-evaluable Mapping. + import numpy as np + from sympde.topology import IdentityMapping, InterfaceMapping + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + assert isinstance(itf.minus, IdentityMapping) + assert isinstance(itf.plus, IdentityMapping) + expected = itf.minus.get_callable_mapping().jacobian(0.3, 0.4) + assert np.array_equal(itf.minus.jacobian(0.3, 0.4), expected) + assert itf.minus(0.3, 0.4) == (0.3, 0.4) + + +def test_copy_preserves_user_set_callable_mapping(): + # post-commit finding A: copy() used to write the name-mangled + # `__callable_map` instead of `_callable_map`, silently dropping a + # user-supplied callable mapping (set via set_callable_mapping) on copy -- + # observable e.g. through InterfaceMapping, which always copies its legs. + # WP13/D2: AnalyticMapping.set_callable_mapping now raises (it is always + # its own callable), so this is retargeted to a plain SymbolicMapping. + class Custom(BasicCallableMapping): + def __call__(self, *eta): return eta + def jacobian(self, *eta): return None + def jacobian_inv(self, *eta): return None + def metric(self, *eta): return None + def metric_det(self, *eta): return None + ldim = 2 + pdim = 2 + + F = SymbolicMapping('F', dim=2) + cm = Custom() + with pytest.warns(DeprecationWarning): + F.set_callable_mapping(cm) + # copy() assigns `_callable_map` directly, not through the deprecated + # setter, so it must not warn. Sitting outside the `pytest.warns` block + # above does NOT assert that -- pytest.warns says nothing about warnings + # raised elsewhere -- so pin it explicitly, as the guards' test does. + with warnings.catch_warnings(): + warnings.simplefilter('error') + G = F.copy() + assert G._callable_map is cm + + +# -- work-package 03 / 06d-4a: InverseMapping / InterfaceMapping / +# MultiPatchMapping are PURE StructuralMapping subclasses (WP03 added +# StructuralMapping alongside Mapping; 06d-4a severed Mapping). + +def test_structural_subclasses_hierarchy(): + from sympde.topology import (InverseMapping, InterfaceMapping, + MultiPatchMapping, Mapping) + for cls in (InverseMapping, InterfaceMapping, MultiPatchMapping): + assert issubclass(cls, StructuralMapping) + assert issubclass(cls, SymbolicMapping) # still a symbolic mapping + assert not issubclass(cls, Mapping) # 06d-4a: Mapping severed + assert not inspect.isabstract(cls) + # own concrete ldim/pdim, not StructuralMapping's abstract ones + assert cls.ldim is not StructuralMapping.ldim + assert cls.pdim is not StructuralMapping.pdim + + +def test_interface_mapping_rejects_point_call_but_stays_domain_callable(): + from sympde.topology import IdentityMapping, InterfaceMapping, Square + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + assert itf.ldim == 2 + with pytest.raises(TypeError, match='StructuralMapping'): + itf(0.3, 0.4) + mapped = itf(Square('D')) + assert type(mapped).__name__ in ('Domain', 'MappedDomain') + + +def test_multi_patch_mapping_rejects_point_call(): + from sympde.topology import MultiPatchMapping, IdentityMapping, Square + mp = MultiPatchMapping({'p1': IdentityMapping('F1', dim=2), + 'p2': IdentityMapping('F2', dim=2)}) + assert mp.ldim == 2 + with pytest.raises(TypeError, match='StructuralMapping'): + mp(0.3, 0.4) + # 06d-4a fixed the parked bug: MultiPatchMapping.__new__ now sets _name, so + # `mp.name`, `mp == mp` and the domain-call path all work. + assert mp.name == 'F1|F2' + assert mp == mp + mapped = mp(Square('D')) + assert type(mapped).__name__ in ('Domain', 'MappedDomain') + + +def test_structural_mapping_still_abstract(): + # WP01 invariant must not regress + assert inspect.isabstract(StructuralMapping) + with pytest.raises(TypeError, match='abstract'): + StructuralMapping('F') + + +# -- post-review fixes (found by /code-review on the WP03 diff) + +def test_structural_mapping_call_accepts_domain_keyword(): + # StructuralMapping.__call__(self, *args) initially dropped the `domain` + # keyword that Mapping.__call__(self, domain) supported -- restored via + # explicit positional/keyword normalization before delegating. + from sympde.topology import IdentityMapping, InterfaceMapping, Square + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + positional = itf(Square('D')) + keyword = itf(domain=Square('D')) + assert type(positional) is type(keyword) + + +def test_structural_mapping_call_returns_mapped_domain(): + # 06d-4a: StructuralMapping.__call__ owns the domain-call body directly + # (Mapping is no longer in the MRO to super() into). Observable result is + # unchanged: a single BasicDomain -> a mapped domain. + from sympde.topology import IdentityMapping, InterfaceMapping, Square + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + assert type(itf(Square('D'))).__name__ in ('Domain', 'MappedDomain') + + +def test_structural_ldim_pdim_are_concrete(): + # 06d-4a: the `ldim = Mapping.ldim` MRO work-around is gone -- each + # structural class has its own concrete ldim/pdim returning stored values. + from sympde.topology import (InverseMapping, InterfaceMapping, + IdentityMapping) + F = IdentityMapping('F', dim=2) + itf = InterfaceMapping(IdentityMapping('A', dim=2), IdentityMapping('B', dim=2)) + assert InverseMapping(F).ldim == 2 and InverseMapping(F).pdim == 2 + assert itf.ldim == 2 and itf.pdim == 2 + assert InterfaceMapping.ldim is not StructuralMapping.ldim + assert InverseMapping.pdim is not StructuralMapping.pdim + + +def test_interface_mapping_copy_reconstructs(): + # 06d-4a fixed the parked bug: InterfaceMapping used to inherit + # Mapping.copy() (type(self)(self.name, ldim=...) -- wrong signature). + # It now has its own copy() that rebuilds from the stored legs. + from sympde.topology import IdentityMapping, InterfaceMapping + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + c = itf.copy() + assert isinstance(c, InterfaceMapping) + assert c == itf + assert c.minus.name == itf.minus.name and c.plus.name == itf.plus.name + + +def test_structural_mappings_are_not_analytic(): + # 06d-4a: the structural object itself carries no analytic `_expressions`. + # 06d-4c: the callable-mapping API lives on SymbolicMapping now (psydac + # attaches spline callables to undefined mappings), so a structural mapping + # inherits get_callable_mapping -- but it has nothing to build / return. + from sympde.topology import IdentityMapping, InterfaceMapping, MultiPatchMapping + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + mp = MultiPatchMapping({'p': IdentityMapping('F', dim=2)}) + for m in (itf, mp): + assert getattr(m, '_expressions', None) is None + with pytest.raises((ValueError, TypeError, AttributeError)): + m.get_callable_mapping() + + +def test_structural_mappings_reject_set_callable_mapping(): + # WP09: a StructuralMapping is symbolic and not point-evaluable (__call__ + # already rejects point calls) -- set_callable_mapping (inherited + # unguarded from SymbolicMapping otherwise) must not be able to silently + # make get_callable_mapping() start returning something. Mirrors + # DiscreteMapping.set_callable_mapping's guard (WP07d). + from sympde.topology import IdentityMapping, InterfaceMapping, MultiPatchMapping + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + mp = MultiPatchMapping({'p': IdentityMapping('F', dim=2)}) + for m in (itf, mp): + with pytest.raises(TypeError): + m.set_callable_mapping(object()) + + +def test_analytic_mapping_set_callable_mapping_raises(): + # WP13/D2: mirrors DiscreteMapping (WP07d) / StructuralMapping (WP09) -- + # an AnalyticMapping is already its own callable mapping. + from sympde.topology import IdentityMapping, PolarMapping + for F in (IdentityMapping('F', dim=2), + PolarMapping('F', dim=2, rmin=0., rmax=1., c1=0., c2=0.)): + with pytest.raises(TypeError, match='to_defined_mapping'): + F.set_callable_mapping(object()) + + +def test_analytic_mapping_get_callable_mapping_is_always_self(): + from sympde.topology import IdentityMapping, PolarMapping + for F in (IdentityMapping('F', dim=2), + PolarMapping('F', dim=2, rmin=0., rmax=1., c1=0., c2=0.)): + assert F.get_callable_mapping() is F + + +def test_symbolic_mapping_set_callable_mapping_is_deprecated_but_works(): + # D3-b: SymbolicMapping.set_callable_mapping is deprecated (conventions: + # deprecated symbols keep working until the explicit remove-aliases + # work-package) -- plain SymbolicMapping.set_callable_mapping/ + # get_callable_mapping stays unguarded, but now warns. + class Custom(BasicCallableMapping): + def __call__(self, *eta): return eta + def jacobian(self, *eta): return None + def jacobian_inv(self, *eta): return None + def metric(self, *eta): return None + def metric_det(self, *eta): return None + ldim = 2 + pdim = 2 + + F = SymbolicMapping('F', dim=2) + cm = Custom() + with pytest.warns(DeprecationWarning): + F.set_callable_mapping(cm) + assert F.get_callable_mapping() is cm + + +def test_set_callable_mapping_warning_names_discrete_mapping(): + # D3-b: the deprecation warning must point callers at the replacement. + class Custom(BasicCallableMapping): + def __call__(self, *eta): return eta + def jacobian(self, *eta): return None + def jacobian_inv(self, *eta): return None + def metric(self, *eta): return None + def metric_det(self, *eta): return None + ldim = 2 + pdim = 2 + + F = SymbolicMapping('F', dim=2) + with pytest.warns(DeprecationWarning, match='DiscreteMapping'): + F.set_callable_mapping(Custom()) + + +def test_set_callable_mapping_guards_do_not_warn(): + # D3-b: the three overriding guards (AnalyticMapping, DiscreteMapping, + # StructuralMapping subclasses) raise TypeError without calling super(), + # so they must not start emitting the new DeprecationWarning. + from sympde.topology import (IdentityMapping, InterfaceMapping, + MultiPatchMapping) + from sympde.topology.mapping import DiscreteMapping + + class Custom(BasicCallableMapping): + def __call__(self, *eta): return eta + def jacobian(self, *eta): return None + def jacobian_inv(self, *eta): return None + def metric(self, *eta): return None + def metric_det(self, *eta): return None + ldim = 2 + pdim = 2 + + F1 = IdentityMapping('F1', dim=2) + F2 = IdentityMapping('F2', dim=2) + G = DiscreteMapping(Custom(), 'G') + itf = InterfaceMapping(F1, F2) + mp = MultiPatchMapping({'p1': F1, 'p2': F2}) + + with warnings.catch_warnings(): + warnings.simplefilter('error') + for mapping in (F1, G, itf, mp): + with pytest.raises(TypeError): + mapping.set_callable_mapping(Custom()) + + +def test_symbolicexpr_lowers_indexed_structural_mapping_to_coordinate(): + # 06d-4a-1: SymbolicExpr.eval's Indexed branch must recognise a severed + # structural mapping (SymbolicMapping, not Mapping) so `itf[i]` lowers to a + # coordinate symbol, not the `{base.name}_{i}` fallback (an invalid + # 'F1|F2_0'). + from sympde.topology import IdentityMapping, InterfaceMapping + from sympde.topology.mapping import SymbolicExpr + from sympy import Symbol + itf = InterfaceMapping(IdentityMapping('F1', dim=2), IdentityMapping('F2', dim=2)) + assert SymbolicExpr(itf[0]) == Symbol('x') + assert SymbolicExpr(itf[1]) == Symbol('y') + + +def test_empty_multipatch_mapping_still_constructs(): + # 06d-4a-1: MultiPatchMapping({}) constructed pre-WP06d-4a via Basic.__new__; + # the _name/_coordinates derivation must not StopIteration on an empty dict. + from sympde.topology import MultiPatchMapping + mp = MultiPatchMapping({}) + assert mp.name == '' + + +def test_mapped_multipatch_domain_with_interface_is_broken_pre_existing(): + # post-review finding: MappedDomain.__new__ rebuilds each patch Interface + # as Interface(e.name, mapping(e.minus), mapping(e.plus)) without + # forwarding e.ornt (or mapping/logical_domain). Since commit 9783335 + # ("Rigorously check arguments of Interface constructor") made + # Interface.__new__ validate ornt strictly, this now raises instead of + # silently building an under-specified Interface. Confirmed pre-existing + # (reproduces identically on the pre-WP03 baseline via git stash) -- + # unrelated to StructuralMapping/Mapping base-list changes, but it lives + # in mapping.py's MappedDomain.__new__, a function this file's tests + # exercise elsewhere. Documenting rather than fixing: the fix belongs to + # whoever owns the Interface/ornt validation, not this work-package. + from sympde.topology import Domain, Square, IdentityMapping + patches = [Square('D1'), Square('D2')] + domain = Domain.join(patches, [((0, 0, 1), (1, 0, -1), 1)], 'domain') + F = IdentityMapping('F', dim=2) + with pytest.raises(AssertionError, match='ornt'): + F(domain) + + +# -- work-package 06a: Mapping re-parented onto SymbolicMapping + +def test_mapping_is_a_symbolic_mapping(): + from sympde.topology import Mapping, IdentityMapping, InterfaceMapping + assert issubclass(Mapping, SymbolicMapping) + # every branch of the hierarchy is now isinstance(_, SymbolicMapping): + with pytest.warns(DeprecationWarning): # 06d-4c + bare = Mapping('F', dim=2) + assert isinstance(bare, SymbolicMapping) # undefined + assert isinstance(IdentityMapping('G', dim=2), SymbolicMapping) # analytic + itf = InterfaceMapping(IdentityMapping('A', dim=2), IdentityMapping('B', dim=2)) + assert isinstance(itf, SymbolicMapping) # structural + # SymbolicMapping appears exactly once in the MRO (no C3 duplication) + assert [c.__name__ for c in Mapping.__mro__].count('SymbolicMapping') == 1 + + +# -- work-package 06c: AnalyticMapping is its own callable mapping + +def test_analytic_mapping_is_its_own_callable_mapping(): + from sympde.topology import IdentityMapping, Square + F = IdentityMapping('F', dim=2) + assert F.get_callable_mapping() is F + # point evaluation (positional) and domain call (positional or `domain=`) + assert F(0.3, 0.4) == (0.3, 0.4) + assert type(F(Square('D'))).__name__ in ('Domain', 'MappedDomain') + assert type(F(domain=Square('D'))).__name__ in ('Domain', 'MappedDomain') + + +def test_analytic_mapping_with_symbolic_constants_rejects_point_call(): + from sympde.topology import PolarMapping + P = PolarMapping('P', dim=2) # rmin/rmax/c1/c2 left symbolic + with pytest.raises(ValueError, match='symbolic constants'): + P(0.5, 0.5) + + +def test_analytic_mapping_rejects_attached_callable(): + # WP13/D2: an AnalyticMapping is already its own callable mapping, so + # attaching a different one (which used to silently "win" for point + # evaluation while is_analytical stayed True for assembly) now raises. + # (Historically -- pre-WP13 -- this test asserted the opposite: that the + # attached callable "won". Kept as its own regression test, alongside the + # more general test_analytic_mapping_set_callable_mapping_raises above, + # because it exercises a real BasicCallableMapping implementation rather + # than a bare object().) + from sympde.topology import PolarMapping + + class Const(BasicCallableMapping): + ldim = pdim = 2 + def __call__(self, *eta): return (7.0, 8.0) + def jacobian(self, *eta): return [[1.0, 0.0], [0.0, 1.0]] + def jacobian_inv(self, *eta): return [[1.0, 0.0], [0.0, 1.0]] + def metric(self, *eta): return [[1.0, 0.0], [0.0, 1.0]] + def metric_det(self, *eta): return 1.0 + + F = PolarMapping('F', dim=2, rmin=0.0, rmax=1.0, c1=0.0, c2=0.0) + with pytest.raises(TypeError, match='to_defined_mapping'): + F.set_callable_mapping(Const()) + + +# -- work-package 06d-2: CallableMapping deleted, BasicMapping folded away + +def test_callable_mapping_is_removed(): + import sympde.topology.callable_mapping as cm + with pytest.raises(AttributeError, match='removed'): + _ = cm.CallableMapping + + +def test_bare_mapping_with_expressions_rejects_get_callable_mapping(): + # a Mapping subclass carrying _expressions but not parented under + # AnalyticMapping is now a mistake -- there is no CallableMapping to build. + from sympde.topology import Mapping + + class M(Mapping): + _expressions = {'x': '2*x1', 'y': '3*x2'} + _ldim = 2 + _pdim = 2 + + with pytest.raises(TypeError, match='AnalyticMapping'): + M('m').get_callable_mapping() + + +def test_undefined_mapping_still_valueerrors_on_get_callable_mapping(): + # 06d-4c: get_callable_mapping() moved onto SymbolicMapping (an undefined + # mapping with no _expressions and no attached callable still ValueErrors). + with pytest.raises(ValueError): + SymbolicMapping('F', dim=2).get_callable_mapping() + + +def test_basicmapping_alias_removed(): + # 06d-2 folded BasicMapping into SymbolicMapping (keeping a deprecated + # alias); 06d-4b drops the alias; 06d-4c relocates SymbolicMapping out of + # sympde.core.basic into sympde.topology.mapping. + with pytest.raises(ImportError): + from sympde.core.basic import BasicMapping as _ # noqa: F401 + with pytest.raises(ImportError): + from sympde.core.basic import SymbolicMapping as _ # noqa: F401 + + from sympde.topology import Mapping, IdentityMapping + assert issubclass(Mapping, SymbolicMapping) + assert isinstance(IdentityMapping('G', dim=2), SymbolicMapping) + # SymbolicMapping appears exactly once in the MRO (BasicMapping is gone) + assert [c.__name__ for c in Mapping.__mro__].count('SymbolicMapping') == 1 + + +# -- work-package 06d-4c: SymbolicMapping is the undefined-mapping constructor, +# `Mapping` is a DeprecationWarning shell over it + +def test_symbolicmapping_is_the_undefined_mapping_constructor(): + from sympde.topology import Square + F = SymbolicMapping('F', dim=2) + assert F.name == 'F' and F.ldim == 2 and F.pdim == 2 + assert F.is_analytical is False + assert type(F(Square('D'))).__name__ in ('Domain', 'MappedDomain') + + +def test_bare_mapping_construction_is_deprecated(): + from sympde.topology import Mapping, IdentityMapping + with pytest.warns(DeprecationWarning, match='SymbolicMapping'): + Mapping('M', dim=2) + with warnings.catch_warnings(): # analytic subclasses do NOT warn + warnings.simplefilter('error', DeprecationWarning) + IdentityMapping('G', dim=2) + + +def test_bare_mapping_rebuild_does_not_warn(): + # /code-review finding 3 (06d-4c-1a): sympy's internal `func(*args)` rebuild + # (Basic.rebuild, cse, pickling, deepcopy) re-enters Mapping.__new__ with + # `name` already a Symbol -- it must stay quiet, or a downstream running + # `filterwarnings=error` breaks on expressions it merely stored / copied. + import copy + import pickle # nosec B403 -- only round-trips objects built by the test + from sympde.topology import Mapping + + with warnings.catch_warnings(): + warnings.simplefilter('ignore', DeprecationWarning) + M = Mapping('F_rebuild', dim=2) + + with warnings.catch_warnings(): + warnings.simplefilter('error', DeprecationWarning) + M.func(*M.args) # the exact rebuild path + copy.deepcopy(M) + pickle.loads(pickle.dumps(M)) # nosec B301 + + +def test_symbolicmapping_honours_injected_jacobian(): + # /code-review finding 2 (06d-4c-1a): the `jacobian=` constructor hook that + # Mapping.__new__ honoured must survive the move into SymbolicMapping.__new__, + # and must not leak into the analytic-constants dict of the Mapping shell. + from sympy import ImmutableDenseMatrix, eye + from sympde.topology import Mapping + + J = ImmutableDenseMatrix(eye(2)) + + F = SymbolicMapping('F_injjac', dim=2, jacobian=J) + assert F.jacobian_symbol is J + + with pytest.warns(DeprecationWarning): # bare Mapping still deprecated + G = Mapping('G_injjac', dim=2, jacobian=J) + # forwarded through the shell to SymbolicMapping.__new__ (not swallowed by + # `**kwargs` and replaced by the default JacobianSymbol). + assert G.jacobian_symbol is J + + +def test_interface_mapping_from_bare_mapping_legs_is_warning_free_on_copy(): + # /code-review finding 2 (06d-4c-1b): SymbolicMapping.copy() re-invokes the + # constructor; for a bare Mapping leg that re-tripped the 06d-4c-1a + # `isinstance(name, str)` warning gate. InterfaceMapping.__new__ copies both + # legs, so building one from bare-Mapping legs must not warn -- while an + # explicit Mapping(...) call still does. + from sympde.topology import Mapping, InterfaceMapping + + with warnings.catch_warnings(): + warnings.simplefilter('ignore', DeprecationWarning) + a = Mapping('a_bare', dim=2) + b = Mapping('b_bare', dim=2) + + with warnings.catch_warnings(): + warnings.simplefilter('error', DeprecationWarning) + InterfaceMapping(a, b) # copies a, b internally -- silent + + with pytest.warns(DeprecationWarning): # explicit construction still warns + Mapping('c_bare', dim=2) + + +# -- WP07: DiscreteMapping -- a DefinedMapping backed by an external callable + +class _FakeCallable(BasicCallableMapping): + """ Minimal BasicCallableMapping for the DiscreteMapping tests. """ + def __init__(self, ldim=2, pdim=2, name=None): + self._l, self._p, self._n = ldim, pdim, name + def __call__(self, *e): return tuple(e) + def jacobian(self, *e): return [[1, 0], [0, 1]] + def jacobian_inv(self, *e): return [[1, 0], [0, 1]] + def metric(self, *e): return [[1, 0], [0, 1]] + def metric_det(self, *e): return 1.0 + @property + def ldim(self): return self._l + @property + def pdim(self): return self._p + @property + def name(self): return self._n + + +def test_discrete_mapping_is_a_concrete_defined_mapping(): + from sympde.topology import DiscreteMapping + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + assert isinstance(G, DefinedMapping) + assert isinstance(G, SymbolicMapping) + assert not inspect.isabstract(DiscreteMapping) + assert G.name == 'D' and G.ldim == 2 and G.pdim == 2 + assert G.is_analytical is False + + +def test_discrete_mapping_delegates_point_evaluation(): + from sympde.topology import DiscreteMapping + + f = _FakeCallable() + G = DiscreteMapping(f, name='D', dim=2) + assert G.get_callable_mapping() is f + assert G.jacobian(0.1, 0.2) == [[1, 0], [0, 1]] + assert G(0.3, 0.4) == (0.3, 0.4) # point call, delegated + + +def test_discrete_mapping_is_callable_on_a_domain(): + from sympde.topology import DiscreteMapping, Square + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + Om = G(Square('S')) + assert type(Om).__name__ in ('Domain', 'MappedDomain') + assert Om.mapping is G + assert Om.logical_domain == Square('S') + + +def test_discrete_mapping_copy_round_trips(): + from sympde.topology import DiscreteMapping + + f = _FakeCallable() + G = DiscreteMapping(f, name='D', dim=2) + c = G.copy() + assert isinstance(c, DiscreteMapping) + assert c == G + assert c.get_callable_mapping() is f + + +def test_discrete_mapping_name_inference_and_error(): + from sympde.topology import DiscreteMapping + + assert DiscreteMapping(_FakeCallable(name='H')).name == 'H' # inferred + with pytest.raises(ValueError): + DiscreteMapping(_FakeCallable()) # no name anywhere + + +def test_discrete_mapping_rejects_non_callable(): + from sympde.topology import DiscreteMapping + with pytest.raises(TypeError): + DiscreteMapping(object(), name='D', dim=2) + + +def test_defined_mapping_factory_sugar(): + from sympde.topology import DiscreteMapping + + G = DefinedMapping(_FakeCallable(), name='D', dim=2) + assert isinstance(G, DiscreteMapping) + # a plain abstract call still fails + with pytest.raises(TypeError, match='abstract'): + DefinedMapping('F') + + +def test_discrete_mapping_usable_as_interface_leg(): + from sympde.topology import DiscreteMapping, InterfaceMapping + + a = DiscreteMapping(_FakeCallable(), name='a', dim=2) + b = DiscreteMapping(_FakeCallable(), name='b', dim=2) + itf = InterfaceMapping(a, b) + assert itf.is_analytical is False + + +# -- WP07-1: post-/code-review fixes for DiscreteMapping + +def test_discrete_mapping_identity_includes_the_wrapped_callable(): + # F1: two same-named DiscreteMappings over different geometries must be + # distinct, else MappedDomain's @cacheit conflates them. + from sympde.topology import DiscreteMapping, Square + + fa, fb = _FakeCallable(), _FakeCallable() + Ga = DiscreteMapping(fa, name='G', dim=2) + Gb = DiscreteMapping(fb, name='G', dim=2) + assert Ga != Gb + assert hash(Ga) != hash(Gb) + assert len({Ga, Gb}) == 2 + + Da, Db = Ga(Square('S1')), Gb(Square('S1')) + assert Da is not Db + assert Da.mapping.get_callable_mapping() is fa + assert Db.mapping.get_callable_mapping() is fb + + # same callable + name + dims -> still equal + assert DiscreteMapping(fa, name='G', dim=2) == Ga + + +def test_discrete_mapping_dims_come_from_the_callable(): + # F4: a dim= / ldim= / pdim= argument may only confirm the callable's dims. + from sympde.topology import DiscreteMapping + + surf = _FakeCallable(ldim=2, pdim=3) + G = DiscreteMapping(surf, name='S', ldim=2, pdim=3) + assert G.ldim == 2 and G.pdim == 3 + + with pytest.raises(ValueError): + DiscreteMapping(surf, name='S2', dim=2) # 2 != pdim 3 + with pytest.raises(ValueError): + DiscreteMapping(_FakeCallable(), name='S3', pdim=5) + + assert DiscreteMapping(_FakeCallable(), name='D', dim=2).pdim == 2 # consistent, ok + + +def test_discrete_mapping_copy_preserves_interface_tags(): + # F2: copy() must carry _is_plus / _is_minus (and coordinates), not just + # (callable, name, ldim, pdim). + from sympde.topology import DiscreteMapping, InterfaceMapping + + itf = InterfaceMapping(DiscreteMapping(_FakeCallable(), 'a', dim=2), + DiscreteMapping(_FakeCallable(), 'b', dim=2)) + assert itf.minus.is_minus is True + assert itf.minus.copy().is_minus is True + assert itf.plus.copy().is_plus is True + + +def test_discrete_mapping_reconstruction_fails_clearly(): + # F3: sympy's func(*args) / pickle / deepcopy shape must raise an actionable + # error, not a confusing "got Symbol". + from sympy import Symbol, Tuple + from sympde.topology import DiscreteMapping + + with pytest.raises(TypeError, match='reconstructed'): + DiscreteMapping(Symbol('D'), Tuple(2)) + + +def test_discrete_mapping_get_callable_guards_none(): + # F3: defensive -- a callable-less DiscreteMapping raises, not AttributeError. + from sympde.topology import DiscreteMapping + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + G._callable_map = None + with pytest.raises(ValueError): + G.get_callable_mapping() + + +def test_discrete_mapping_has_callable_mapping(): + # WP07e: a non-raising predicate for get_callable_mapping()'s guard, so + # callers don't need to catch ValueError as control flow. + from sympde.topology import DiscreteMapping + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + assert G.has_callable_mapping() is True + + G._callable_map = None + assert G.has_callable_mapping() is False + + +def test_discrete_mapping_set_callable_mapping_raises(): + # WP07d: the wrapped callable is fixed at construction (part of identity, + # see _hashable_content) -- set_callable_mapping must not be able to swap + # it, unlike the base SymbolicMapping.set_callable_mapping it would + # otherwise inherit unguarded. + from sympde.topology import DiscreteMapping + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + with pytest.raises(TypeError): + G.set_callable_mapping(_FakeCallable()) + + +# -- WP07-2: second /code-review round for DiscreteMapping + +def test_discrete_mapping_jacobian_expr_uses_the_final_identity(): + # G1: __new__ must attach _callable_map before building _jac / _metric, so + # `G[i]` and the Indexed(G, i) baked into jacobian_expr hash identically. + from sympy import Symbol + from sympde.topology import DiscreteMapping + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + assert G[0] in G.jacobian_expr.free_symbols + assert G[1] in G.metric_expr.free_symbols + assert G.jacobian_expr.subs(G[0], Symbol('v')) != G.jacobian_expr + + Ga = DiscreteMapping(_FakeCallable(), name='S', dim=2) + Gb = DiscreteMapping(_FakeCallable(), name='S', dim=2) + assert all(s.base is Ga for s in Ga.jacobian_expr.free_symbols + if hasattr(s, 'base')) + assert all(s.base is Gb for s in Gb.jacobian_expr.free_symbols + if hasattr(s, 'base')) + + +def test_discrete_mapping_delegators_guard_none_callable(): + # G2: every delegating point-eval method routes through + # get_callable_mapping(), so a None _callable_map gives a ValueError, not + # AttributeError / TypeError. + from sympde.topology import DiscreteMapping + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + G._callable_map = None + for call in (lambda: G.jacobian(0., 0.), + lambda: G.jacobian_inv(0., 0.), + lambda: G.metric(0., 0.), + lambda: G.metric_det(0., 0.), + lambda: G(0.1, 0.2)): + with pytest.raises(ValueError): + call() + + +def test_discrete_mapping_name_rules(): + # G3: an explicit name wins and must be non-empty; a falsy name is not + # silently replaced by the callable's. + from sympde.topology import DiscreteMapping + + with pytest.raises(ValueError): + DiscreteMapping(_FakeCallable(name='H'), name='') + with pytest.raises(ValueError): + DiscreteMapping(_FakeCallable(name=None)) + assert DiscreteMapping(_FakeCallable(name='H')).name == 'H' + assert DiscreteMapping(_FakeCallable(name='H'), name='K').name == 'K' + + +def test_discrete_mapping_rejects_unexpected_kwargs(): + from sympde.topology import DiscreteMapping + with pytest.raises(TypeError): + DiscreteMapping(_FakeCallable(), name='D', dim=2, bogus=1) + + +# -- WP07b-1: DiscreteMapping survives sympy's func(*args) rebuild + +def test_discrete_mapping_func_is_the_identity_rebuild(): + # sympy walks an expression tree calling `node.func(*node.args)` + # (Basic.rebuild, cse, and sympy.simplify's replace-reducer). The base + # `Basic.func` is `type(self)`, so it would call + # `DiscreteMapping(Symbol(name), Tuple(pdim))` -> TypeError. A DiscreteMapping + # is a symbolic leaf (the callable is not in `.args`), so `func(*args)` must + # be the identity -- this is what lets a bilinear form be discretised on a + # DiscreteMapping-carried domain (Jacobian(M).inv() -> simplify -> func). + import sympy + from sympde.topology import DiscreteMapping + from sympde.topology.mapping import Jacobian + + G = DiscreteMapping(_FakeCallable(), name='D', dim=2) + + assert G.func(*G.args) is G + assert G.func(sympy.Symbol('x'), (99,), foo='bar') is G + assert sympy.simplify(G[0]**2 / G[1]) == G[0]**2 / G[1] # no crash + Jacobian(G).inv() # the assembly path; no crash diff --git a/sympde/topology/tests/test_space.py b/sympde/topology/tests/test_space.py index 793afe55..6f16d16c 100644 --- a/sympde/topology/tests/test_space.py +++ b/sympde/topology/tests/test_space.py @@ -9,6 +9,35 @@ from sympde.topology import ScalarFunction, VectorFunction from sympde.topology import Projector + +def test_function_symbol_identity(): + domain = Domain('Omega', dim=2) + scalar_space = ScalarFunctionSpace('V', domain) + equivalent_scalar_space = ScalarFunctionSpace('V', domain) + other_scalar_space = ScalarFunctionSpace('W', domain) + vector_space = VectorFunctionSpace('X', domain) + other_vector_space = VectorFunctionSpace('Y', domain) + + u = ScalarFunction(scalar_space, 'u') + same_u = ScalarFunction(equivalent_scalar_space, 'u') + other_u = ScalarFunction(other_scalar_space, 'u') + v = VectorFunction(vector_space, 'v') + other_v = VectorFunction(other_vector_space, 'v') + + assert scalar_space == equivalent_scalar_space + assert hash(scalar_space) == hash(equivalent_scalar_space) + assert scalar_space != other_scalar_space + assert u == same_u + assert hash(u) == hash(same_u) + assert u != other_u + assert v != other_v + assert u + other_u != 2 * u + assert (u + other_u).free_symbols == {u, other_u} + assert u.is_commutative + assert not v.is_commutative + assert u.free_symbols == {u} + assert v.free_symbols == {v} + #============================================================================== def test_space_1d_1(): diff --git a/sympde/topology/tests/test_symbolic_expr.py b/sympde/topology/tests/test_symbolic_expr.py index 43d1151c..b49e14f6 100644 --- a/sympde/topology/tests/test_symbolic_expr.py +++ b/sympde/topology/tests/test_symbolic_expr.py @@ -4,7 +4,7 @@ from sympde.topology import Domain, ScalarFunctionSpace, element_of from sympde.topology import dx, dy, dz from sympde.topology import dx1, dx2, dx3 -from sympde.topology import Mapping, Jacobian +from sympde.topology import SymbolicMapping, Jacobian from sympde.topology import LogicalExpr from sympde.topology import SymbolicExpr @@ -43,7 +43,7 @@ def test_derivatives_2d_without_mapping(): def test_derivatives_2d_with_mapping(): dim = 2 - M = Mapping('M', dim=dim) + M = SymbolicMapping('M', dim=dim) O = M(Domain('Omega', dim=dim)) V = ScalarFunctionSpace('V', O, kind='h1') u = element_of(V, 'u') diff --git a/sympde/topology/tests/test_topology.py b/sympde/topology/tests/test_topology.py index 704affc0..892edc69 100644 --- a/sympde/topology/tests/test_topology.py +++ b/sympde/topology/tests/test_topology.py @@ -1,9 +1,12 @@ import os +from sympy import preorder_traversal +from sympy.core import Basic + from sympde.topology import InteriorDomain, Union from sympde.topology import Boundary from sympde.topology import Domain, ElementDomain -from sympde.topology import Area, Mapping +from sympde.topology import Area, SymbolicMapping from sympde.topology import Interface from sympde.topology import Line, Square, Cube from sympde.topology import IdentityMapping @@ -11,6 +14,33 @@ base_dir = os.path.dirname(os.path.realpath(__file__)) topo_dir = os.path.join(base_dir, 'data') + +def test_structural_arguments_are_sympy_objects(): + domain = Domain('Omega', dim=2) + boundary = Boundary(r'\Gamma_1', domain) + interface = Interface('I', boundary, boundary, ornt=1) + square = Square('S', bounds1=(-1, 2), bounds2=(3, 4)) + + assert domain.name == 'Omega' + assert boundary.name == r'\Gamma_1' + assert boundary.axis is None + assert interface.ornt == 1 + + for expr in (domain, domain.interior, boundary, interface, square, + square.interior, square.boundary.args[0]): + assert all(isinstance(arg, Basic) for arg in preorder_traversal(expr)) + assert expr.func(*expr.args) == expr + + replacement = Domain('Replacement', dim=2) + replaced_boundary = boundary.xreplace({domain: replacement}) + assert replaced_boundary.domain == replacement + + replacement_boundary = Boundary(r'\Gamma_2', domain) + replaced_interface = interface.xreplace({boundary: replacement_boundary}) + assert replaced_interface.minus == replacement_boundary + assert replaced_interface.plus == replacement_boundary + + #============================================================================== def test_interior_domain(): D1 = InteriorDomain('D1', dim=2) @@ -36,8 +66,8 @@ def test_topology_1(): B = Square('B') # ... - M1 = Mapping('M1', dim=2) - M2 = Mapping('M2', dim=2) + M1 = SymbolicMapping('M1', dim=2) + M2 = SymbolicMapping('M2', dim=2) D1 = M1(A) D2 = M2(B)