New mapping class hierarchy (sympde part)

[campospinto]

Summary

This PR replaces the single Mapping class with a small hierarchy that keeps
two questions apart: can this mapping be used symbolically? and can it be evaluated at a point?
Today Mapping answers both. It is symbolic, and it
becomes point-evaluable only when it has _expressions or when a numerical
mapping is attached to it after construction with set_callable_mapping.

The new root is SymbolicMapping. Under it:

• DefinedMapping is for mappings that are symbolic and point-evaluable,
  with two concrete classes:
  • AnalyticMapping for mappings given by formulas,
  • DiscreteMapping for mappings backed by a numerical mapping, such as a
    psydac spline.
• StructuralMapping is for mappings that are symbolic only: InverseMapping,
  InterfaceMapping and MultiPatchMapping.

The main user-facing gain: a spline geometry can now have its own symbolic
identity (DiscreteMapping(F_h, 'F')). You can build a domain from it, join
patches, and discretize, without mutating another object.

Companion PR: pyccel/psydac#TBD. It uses the new classes and depends on this
PR. sympde itself has no dependency on psydac.

Main changes

New classes (sympde/topology/mapping.py, exported from sympde.topology)

• SymbolicMapping(IndexedBase): the new root. SymbolicMapping('F', dim=2)
  is an undefined mapping. You can apply it to a domain and use it in forms,
  and it carries the symbolic Jacobian and metric. It cannot be evaluated at a
  point.
• DefinedMapping(SymbolicMapping, BasicCallableMapping): an abstract class for
  mappings that are symbolic and point-evaluable. As a shortcut,
  DefinedMapping(callable, name) returns a DiscreteMapping.
• AnalyticMapping(DefinedMapping): the base class for mappings given by
  formulas. It turns its own _expressions into NumPy functions the first time
  they are needed, and evaluates both scalars and arrays. It is its own
  numerical mapping, so get_callable_mapping() returns self. All 11 gallery
  mappings now derive from it (IdentityMapping, PolarMapping,
  TargetMapping, CzarnyMapping, TorusMapping, …).
• DiscreteMapping(DefinedMapping): gives a symbolic name to any external
  BasicCallableMapping, such as psydac's SplineCallableMapping, and
  forwards every point evaluation to it. is_analytical is False, so
  assembly evaluates the geometry on the quadrature grid. The wrapped mapping
  is fixed when the object is built and is part of its identity (it counts in
  == and hash). has_callable_mapping() is a check that never raises.
• StructuralMapping(SymbolicMapping): an abstract class for the multipatch
  relations InverseMapping, InterfaceMapping and MultiPatchMapping. You
  can apply them to a domain, but evaluating them at a point raises
  TypeError. They also gain a working .name and .copy().
• BasicCallableMapping (the plain-Python interface for point evaluation) is
  unchanged.

Removed or deprecated

• CallableMapping is removed, because an AnalyticMapping is its own
  numerical mapping. sympde.topology.callable_mapping now only re-exports
  BasicCallableMapping.
• sympde.core.basic.BasicMapping is removed and merged into SymbolicMapping.
• Mapping is kept for now, but constructing it directly emits a
  DeprecationWarning. Legacy subclasses of the form class X(Mapping) with
  _expressions keep working.
• set_callable_mapping is deprecated on SymbolicMapping and raises
  TypeError on AnalyticMapping, DiscreteMapping and StructuralMapping.
• The symbolic property Mapping.jacobian is deprecated in favour of
  jacobian_symbol.

Behaviour changes to watch for downstream

• isinstance(x, Mapping) is now False for analytic and structural mappings.
  Use isinstance(x, SymbolicMapping) to mean "any symbolic mapping".
• On an AnalyticMapping, F.jacobian is now the numerical method
  F.jacobian(*eta). The symbolic Jacobian is F.jacobian_symbol.
• Creating an AnalyticMapping subclass that has no _expressions raises
  TypeError. Use SymbolicMapping(name, dim=...) for an undefined mapping
  instead.
• Domain.from_file builds its patches from SymbolicMapping objects.
  Interface accepts any SymbolicMapping.

Bug fix

• sympde/topology/derivatives.py: the chain rule for logical derivatives used
  only one of the mappings in an expression, picked arbitrarily from a set. It
  now adds up the contributions of all of them.

Tests

• The new sympde/topology/tests/test_mapping_hierarchy.py (about 1,080 lines)
  covers the hierarchy, point evaluation, identity and hashing, the guards and
  the deprecations. Existing tests are migrated to the new names.
• Locally, pytest -n auto --dist loadgroup --pyargs sympde gives
  251 passed, 2 skipped. The examples below run without warnings under
  -W error::DeprecationWarning.

Documentation

Before: the mapping classes on master (sympde) and devel (psydac)

classDiagram
    class IndexedBase {
      <<sympy>>
    }
    class BasicMapping {
      <<sympde.core, empty>>
    }
    class Mapping {
      name ldim pdim _expressions
      __call__(domain) MappedDomain
      jacobian  jacobian_expr  metric_expr
      get_callable_mapping()
      set_callable_mapping(F)
    }
    class InverseMapping
    class InterfaceMapping
    class MultiPatchMapping
    class PolarMapping {
      <<analytic gallery>>
    }
    class BasicCallableMapping {
      <<abc.ABC>>
      __call__(*eta)
      jacobian(*eta)  jacobian_inv(*eta)
      metric(*eta)  metric_det(*eta)
    }
    class CallableMapping {
      lambdified _expressions
    }
    class SplineMapping {
      <<psydac>>
    }
    class NurbsMapping {
      <<psydac>>
    }

    IndexedBase <|-- BasicMapping
    BasicMapping <|-- Mapping
    Mapping <|-- InverseMapping
    Mapping <|-- InterfaceMapping
    Mapping <|-- MultiPatchMapping
    Mapping <|-- PolarMapping
    BasicCallableMapping <|-- CallableMapping
    BasicCallableMapping <|-- SplineMapping
    SplineMapping <|-- NurbsMapping
    Mapping ..> CallableMapping : builds from _expressions
    Mapping o-- BasicCallableMapping : set_callable_mapping

Loading

The problems with this design:

• One class, three roles. Mapping stands for a formula, a placeholder for
  a mapping defined elsewhere, and a relation between patches. Nothing in the
  type says which one you have. InverseMapping, InterfaceMapping and
  MultiPatchMapping inherit a point-evaluation API they cannot honour.
• Point evaluation depends on hidden state. A Mapping without
  _expressions becomes point-evaluable only after something calls
  set_callable_mapping on it. psydac's Geometry.read() did this to the
  mappings returned by Domain.from_file, so whether a mapping could be
  evaluated depended on which code had run before.
• Splines have no symbolic identity. psydac's SplineMapping is purely
  numerical. To use one in a symbolic domain, it had to be attached to some
  other Mapping.

After: the new hierarchy

classDiagram
    class IndexedBase {
      <<sympy>>
    }
    class SymbolicMapping {
      name  ldim  pdim  coordinates  logical_coordinates
      __call__(domain) MappedDomain
      jacobian_symbol  jacobian_expr  metric_expr  is_analytical
      get_callable_mapping()  copy()
    }
    class BasicCallableMapping {
      <<abc.ABC>>
      __call__(*eta)  jacobian(*eta)  jacobian_inv(*eta)
      metric(*eta)  metric_det(*eta)  ldim  pdim
    }
    class DefinedMapping {
      <<abstract>>
      symbolic + point-evaluable
    }
    class AnalyticMapping {
      _expressions, lambdified on demand
    }
    class DiscreteMapping {
      wraps a BasicCallableMapping
      has_callable_mapping()
    }
    class StructuralMapping {
      <<abstract>>
      symbolic only, no point evaluation
    }
    class InverseMapping
    class InterfaceMapping
    class MultiPatchMapping
    class PolarMapping {
      <<analytic gallery>>
    }
    class Mapping {
      <<deprecated>>
    }
    class SplineCallableMapping {
      <<psydac>>
      from_mapping(F, space=None, ...)
      to_defined_mapping(name)
    }
    class NurbsCallableMapping {
      <<psydac>>
    }

    IndexedBase <|-- SymbolicMapping
    SymbolicMapping <|-- DefinedMapping
    BasicCallableMapping <|-- DefinedMapping
    SymbolicMapping <|-- StructuralMapping
    SymbolicMapping <|-- Mapping
    DefinedMapping <|-- AnalyticMapping
    DefinedMapping <|-- DiscreteMapping
    AnalyticMapping <|-- PolarMapping
    StructuralMapping <|-- InverseMapping
    StructuralMapping <|-- InterfaceMapping
    StructuralMapping <|-- MultiPatchMapping
    BasicCallableMapping <|-- SplineCallableMapping
    SplineCallableMapping <|-- NurbsCallableMapping
    DiscreteMapping o-- BasicCallableMapping : wraps

Loading

The psydac classes are shown for context; they come in the companion PR
(SplineMapping is renamed SplineCallableMapping there, and the old name
still works).

class	what it is	can apply to a domain	can evaluate at a point
SymbolicMapping	undefined mapping: a name and its dimensions	yes	no
AnalyticMapping (+ gallery)	mapping given by formulas (_expressions)	yes	yes
DiscreteMapping	symbolic name for an external BasicCallableMapping	yes	yes (forwarded)
InverseMapping, InterfaceMapping, MultiPatchMapping	relations between patches, built by Domain.join	yes	no (TypeError)
BasicCallableMapping	numerical interface only (e.g. psydac splines)	no	yes

Example

This runs with sympde alone. Outputs are shown as comments.

import numpy as np
from sympde.topology import (Square, Domain, SymbolicMapping, AnalyticMapping,
                             DiscreteMapping, BasicCallableMapping, PolarMapping)

A = Square('A', bounds1=(0., 1.), bounds2=(0., 0.5 * np.pi))

# 1. Undefined mapping: symbolic only
F = SymbolicMapping('F', dim=2)
F(A)                        # mapped Domain, usable in forms; .mapping is F
F.jacobian_symbol           # Jacobian(F)
# F.get_callable_mapping()  # ValueError: no point evaluation available

# 2. Analytic mapping: symbolic + point-evaluable (it is its own callable)
P = PolarMapping('P', dim=2, c1=0., c2=0., rmin=0.3, rmax=1.0)
P(A)                        # mapped Domain
P(0.5, 0.7)                 # (0.4971..., 0.4187...)
P.metric_det(0.5, 0.7)      # 0.207025
P.jacobian_expr             # Matrix([[0.7*cos(x2), -(0.7*x1 + 0.3)*sin(x2)], ...])
E1, E2 = np.meshgrid(np.linspace(0, 1, 3), np.linspace(0, 1, 4), indexing='ij')
X, Y = P(E1, E2)            # arrays work too: X.shape == (3, 4)
P.get_callable_mapping() is P   # True

# A new analytic mapping only needs its expressions (constants become kwargs)
class ShearMapping(AnalyticMapping):
    _expressions = {'x': 'x1 + s*x2', 'y': 'x2'}
    _ldim = 2
    _pdim = 2

S = ShearMapping('S', s=0.5)
S(1.0, 2.0)                 # (2.0, 2.0)

# 3. Discrete mapping: give a symbolic identity to an external numeric mapping
class ShearCallable(BasicCallableMapping):
    """Plain numeric mapping (in practice: a psydac SplineCallableMapping)."""
    def __init__(self, s):   self._s = s
    def __call__(self, *eta):     return (eta[0] + self._s * eta[1], eta[1])
    def jacobian(self, *eta):     return np.array([[1., self._s], [0., 1.]])
    def jacobian_inv(self, *eta): return np.linalg.inv(self.jacobian(*eta))
    def metric(self, *eta):       J = self.jacobian(*eta); return J.T @ J
    def metric_det(self, *eta):   return np.linalg.det(self.metric(*eta))
    ldim = property(lambda self: 2)
    pdim = property(lambda self: 2)

G_h = ShearCallable(0.5)
G = DiscreteMapping(G_h, 'G')   # also: DefinedMapping(G_h, 'G')
G(A)                        # mapped Domain; .mapping is G
G(1.0, 2.0)                 # (2.0, 2.0), delegated to G_h
G.is_analytical             # False: assembly evaluates G_h on the grid
G.get_callable_mapping() is G_h   # True
# G.set_callable_mapping(...)     # TypeError: the callable is part of G's identity

# 4. Structural mappings: built by Domain.join, symbolic only
A1 = Square('A1', bounds1=(0., 1.), bounds2=(0., 0.5 * np.pi))
A2 = Square('A2', bounds1=(0., 1.), bounds2=(np.pi, 1.5 * np.pi))
P1 = PolarMapping('P1', dim=2, c1=0.,  c2=0., rmin=0.3, rmax=1.0)
P2 = PolarMapping('P2', dim=2, c1=1.3, c2=0., rmin=1.0, rmax=0.3)
Omega = Domain.join([P1(A1), P2(A2)], [((0, 1, -1), (1, 1, -1), 1)], 'Omega')
type(Omega.mapping).__name__             # 'MultiPatchMapping'
type(Omega.interfaces.mapping).__name__  # 'InterfaceMapping'
# Omega.mapping(0.5, 0.5)                # TypeError: not point-evaluable

With psydac (companion PR), G_h is a spline and step 3 reads
SplineCallableMapping.from_mapping(P, ncells=..., degree=...).to_defined_mapping('G').
discretize(domain) then builds the discrete geometry straight from the
splines, so no geometry file and no ncells are needed.

Migration

before	after
Mapping('F', dim=2)	SymbolicMapping('F', dim=2)
class M(Mapping): _expressions = {...}	class M(AnalyticMapping): _expressions = {...}
isinstance(x, Mapping)	isinstance(x, SymbolicMapping)
F.jacobian (symbolic)	F.jacobian_symbol
F.get_callable_mapping() on an analytic mapping	still works and returns F; F(*eta), F.jacobian(*eta) can be called directly
CallableMapping(F)	F itself (an AnalyticMapping)
M.set_callable_mapping(F_h)	DiscreteMapping(F_h, 'M')
from sympde.core.basic import BasicMapping	from sympde.topology import SymbolicMapping

Mapping, SymbolicMapping.set_callable_mapping and the deprecated
Mapping.jacobian are kept, with warnings, so downstream code has time to
migrate. They will be removed in a follow-up PR.