Add BasicCallableMapping abstract base class (#119)

Add new class to module `sympde.topology.mapping`:
   - Require methods `__call__`, `jacobian`, `metric`, `metric_det`, all called with arguments `(*eta)`
   - Require properties `ldim` and `pdim`
   - Require method `jacobian_inv`, which should raise exception if Jacobian matrix is not invertible (e.g. when `ldim < pdim`)
   - Make `CallableMapping` subclass it
   - Spline and NURBS mappings in Psydac will subclass it

Add new method `set_callable_mapping` to symbolic `Mapping` class, in order to attach any object of type `BasicCallableMapping` to it. Psydac will use it to attach spline and NURBS mappings to MappedDomains with undefined Mapping

Further changes:
   - Allow for mappings with `pdim > ldim` (fixes #70)
   - Update library version to 0.16.0

sympde/topology/analytical_mapping.py

@@ -71,8 +71,6 @@ class CollelaMapping2D(Mapping):
     """
     Represents a Collela 2D Mapping object.
 
-    Examples
-
     """
     _expressions = {'x': '2.*(x1 + eps*sin(2.*pi*k1*x1)*sin(2.*pi*k2*x2)) - 1.',
                     'y': '2.*(x2 + eps*sin(2.*pi*k1*x1)*sin(2.*pi*k2*x2)) - 1.'}
@@ -83,32 +81,65 @@ class CollelaMapping2D(Mapping):
 #==============================================================================
 class TorusMapping(Mapping):
     """
-    Represents a Torus 3D Mapping object.
+    Parametrization of a torus (or a portion of it) of major radius R0, using
+    toroidal coordinates (x1, x2, x3) = (r, theta, phi), where:
 
-    Examples
+      - minor radius    0 <= r < R0
+      - poloidal angle  0 <= theta < 2 pi
+      - toroidal angle  0 <= phi < 2 pi
 
     """
-    _expressions = {'x': '(R0+x1*cos(x2))*cos(x3)',
-                    'y': '(R0+x1*cos(x2))*sin(x3)',
-                    'z': 'x1*sin(x2)'}
+    _expressions = {'x': '(R0 + x1 * cos(x2)) * cos(x3)',
+                    'y': '(R0 + x1 * cos(x2)) * sin(x3)',
+                    'z':       'x1 * sin(x2)'}
 
-    _ldim        = 2
-    _pdim        = 2
+    _ldim        = 3
+    _pdim        = 3
 
 #==============================================================================
-class TwistedTargetMapping(Mapping):
+# TODO [YG, 07.10.2022]: add test in sympde/topology/tests/test_logical_expr.py
+class TorusSurfaceMapping(Mapping):
     """
-    Represents a Twisted Target 3D Mapping object.
+    3D surface obtained by "slicing" the torus above at r = a.
+    The parametrization uses the coordinates (x1, x2) = (theta, phi), where:
 
-    Examples
+      - poloidal angle  0 <= theta < 2 pi
+      - toroidal angle  0 <= phi < 2 pi
 
     """
-    _expressions = {'x': 'c1 + (1-k)*x1*cos(x2) - D*x1**2',
-                    'y': 'c2 + (1+k)*x1*sin(x2)',
-                    'z': 'c3 + x3*x1**2*sin(2*x2)'}
+    _expressions = {'x': '(R0 + a * cos(x1)) * cos(x2)',
+                    'y': '(R0 + a * cos(x1)) * sin(x2)',
+                    'z':       'a * sin(x1)'}
 
     _ldim        = 2
-    _pdim        = 2
+    _pdim        = 3
+
+#==============================================================================
+# TODO [YG, 07.10.2022]: add test in sympde/topology/tests/test_logical_expr.py
+class TwistedTargetSurfaceMapping(Mapping):
+    """
+    3D surface obtained by "twisting" the TargetMapping out of the (x, y) plane
+
+    """
+    _expressions = {'x': 'c1 + (1-k) * x1 * cos(x2) - D *x1**2',
+                    'y': 'c2 + (1+k) * x1 * sin(x2)',
+                    'z': 'c3 + x1**2 * sin(2*x2)'}
+
+    _ldim        = 2
+    _pdim        = 3
+
+#==============================================================================
+class TwistedTargetMapping(Mapping):
+    """
+    3D volume obtained by "extruding" the TwistedTargetSurfaceMapping along z.
+
+    """
+    _expressions = {'x': 'c1 + (1-k) * x1 * cos(x2) - D * x1**2',
+                    'y': 'c2 + (1+k) * x1 * sin(x2)',
+                    'z': 'c3 + x3 * x1**2 * sin(2*x2)'}
+
+    _ldim        = 3
+    _pdim        = 3
 
 #==============================================================================
 class SphericalMapping(Mapping):

sympde/topology/callable_mapping.py

@@ -1,8 +1,12 @@
-from sympde.utilities.utils import lambdify_sympde
-from .mapping import Mapping
 from sympy import Symbol
 
-class CallableMapping:
+from sympde.utilities.utils import lambdify_sympde
+from .mapping import Mapping, BasicCallableMapping
+
+__all__ = ('CallableMapping',)
+
+#==============================================================================
+class CallableMapping(BasicCallableMapping):
 
     def __init__( self, mapping, **kwargs ):
 
@@ -54,27 +58,28 @@ def __init__( self, mapping, **kwargs ):
     #--------------------------------------------------------------------------
     # Abstract interface
     #--------------------------------------------------------------------------
-    def __call__( self, *eta ):
+    def __call__(self, *eta):
         return tuple( f( *eta ) for f in self._func_eval)
 
-    def jacobian( self, *eta ):
+    def jacobian(self, *eta):
         return self._jacobian( *eta )
 
-    def jacobian_inv( self, *eta ):
-        return self._jacobian_inv( *eta )
+    def jacobian_inv(self, *eta):
+        """ Compute the inverse Jacobian matrix, if possible."""
+        return self._jacobian_inv(*eta)
 
-    def metric( self, *eta ):
+    def metric(self, *eta):
         return self._metric( *eta )
 
-    def metric_det( self, *eta ):
+    def metric_det(self, *eta):
         return self._metric_det( *eta )
 
     @property
-    def ldim( self ):
+    def ldim(self):
         return self.symbolic_mapping.ldim
 
     @property
-    def pdim( self ):
+    def pdim(self):
         return self.symbolic_mapping.pdim
 
     #--------------------------------------------------------------------------

sympde/topology/mapping.py

@@ -1,7 +1,5 @@
 # coding: utf-8
-
-
-
+from abc import ABC, abstractmethod
 from sympy                 import Indexed, IndexedBase, Idx
 from sympy                 import Matrix, ImmutableDenseMatrix
 from sympy                 import Function, Expr
@@ -10,7 +8,6 @@
 from sympy.core            import Basic
 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                 import S
@@ -26,7 +23,6 @@
 from sympde.calculus.core     import grad, div, curl, laplace #, hessian
 from sympde.calculus.core     import dot, inner, outer, _diff_ops
 from sympde.calculus.core     import has, DiffOperator
-
 from sympde.calculus.matrices import MatrixSymbolicExpr, MatrixElement, SymbolicTrace, Inverse
 from sympde.calculus.matrices import SymbolicDeterminant, Transpose
 
@@ -48,6 +44,7 @@
 # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping
 
 __all__ = (
+    'BasicCallableMapping',
     'Contravariant',
     'Covariant',
     'InterfaceMapping',
@@ -87,6 +84,52 @@ def get_logical_test_function(u):
     el              = l_space.element(u.name)
     return el
 
+#==============================================================================
+class BasicCallableMapping(ABC):
+    """
+    Transformation of coordinates, which can be evaluated.
+
+    F: R^l -> R^p
+    F(eta) = x
+
+    with l <= p
+    """
+    @abstractmethod
+    def __call__(self, *eta):
+        """ Evaluate mapping at location eta. """
+
+    @abstractmethod
+    def jacobian(self, *eta):
+        """ Compute Jacobian matrix at location eta. """
+
+    @abstractmethod
+    def jacobian_inv(self, *eta):
+        """ Compute inverse Jacobian matrix at location eta.
+            An exception should be raised if the matrix is singular.
+        """
+
+    @abstractmethod
+    def metric(self, *eta):
+        """ Compute components of metric tensor at location eta. """
+
+    @abstractmethod
+    def metric_det(self, *eta):
+        """ Compute determinant of metric tensor at location eta. """
+
+    @property
+    @abstractmethod
+    def ldim(self):
+        """ Number of logical/parametric dimensions in mapping
+            (= number of eta components).
+        """
+
+    @property
+    @abstractmethod
+    def pdim(self):
+        """ Number of physical dimensions in mapping
+            (= number of x components).
+        """
+
 #==============================================================================
 class Mapping(BasicMapping):
     """
@@ -192,10 +235,10 @@ def __new__(cls, name, dim=None, **kwargs):
 
             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()
+                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()
+                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)
@@ -212,6 +255,31 @@ def __new__(cls, name, dim=None, **kwargs):
 
         return obj
 
+    #--------------------------------------------------------------------------
+    # Callable mapping
+    #--------------------------------------------------------------------------
+    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
+
+    def set_callable_mapping(self, F):
+
+        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 ):
         return self._name
@@ -324,13 +392,6 @@ def _hashable_content(self):
                 self._expressions, self._constants, self._is_plus, self._is_minus)
         return tuple([a for a in args if a is not None])
 
-    def get_callable_mapping( self ):
-
-        if self._callable_map is None:
-            import sympde.topology.callable_mapping as cm
-            self._callable_map = cm.CallableMapping( self )
-        return self._callable_map
-
     def _eval_subs(self, old, new):
         return self
 
@@ -1366,4 +1427,3 @@ def eval(cls, *_args, **kwargs):
         # TODO: check if we should use 'sympy.sympify(expr)' instead
         else:
             raise NotImplementedError('Cannot translate to Sympy: {}'.format(expr))
-

sympde/topology/tests/test_callable_mapping.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+from sympde.topology.mapping            import Mapping, BasicCallableMapping
 from sympde.topology.analytical_mapping import IdentityMapping, AffineMapping
 from sympde.topology.analytical_mapping import PolarMapping
 
@@ -466,3 +467,54 @@ def test_polar_mapping_array():
     assert np.allclose(f.jacobian_inv(xx1, xx2), J_inv, rtol=RTOL, atol=ATOL)
     assert np.allclose(f.metric      (xx1, xx2), G    , rtol=RTOL, atol=ATOL)
     assert np.allclose(f.metric_det  (xx1, xx2), G_det, rtol=RTOL, atol=ATOL)
+
+#==============================================================================
+def test_user_defined_callable_mapping():
+
+    class UserIdentity(BasicCallableMapping):
+        """ Identity in N dimensions.
+        """
+
+        def __init__(self, ndim):
+            self._ndim = ndim
+
+        def __call__(self, *eta):
+            assert len(eta) == self._ndim
+            return eta
+
+        def jacobian(self, *eta):
+            assert len(eta) == self._ndim
+            return np.eye(self._ndim)
+
+        def jacobian_inv(self, *eta):
+            assert len(eta) == self._ndim
+            return np.eye(self._ndim)
+
+        def metric(self, *eta):
+            assert len(eta) == self._ndim
+            return np.eye(self._ndim)
+
+        def metric_det(self, *eta):
+            assert len(eta) == self._ndim
+            return 1.0
+
+        @property
+        def ldim(self):
+            return self._ndim
+
+        @property
+        def pdim(self):
+            return self._ndim
+
+    F = Mapping('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
+
+    assert F.get_callable_mapping() is f
+    assert f(4, 5, 6) == (4, 5, 6)
+    assert np.array_equal(f.jacobian(1, 2, 3)     , [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    assert np.array_equal(f.jacobian_inv(-7, 0, 7), [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    assert np.array_equal(f.metric(0, 1, 1)       , [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    assert f.metric_det(-5, 8, -9) == 1.0
+    assert f.ldim == 3
+    assert f.pdim == 3

sympde/version.py

@@ -1 +1 @@
-__version__ = "0.15.2"
+__version__ = "0.16.0"