diff --git a/.gitignore b/.gitignore
index 27146ff3..4227e33b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -111,3 +111,6 @@ doc/utils/*.png
 
 # Swap files
 *.swp
+
+__gpyccel__/*
+__psydac__/*
diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py
index e75189be..7c3c2ee7 100644
--- a/sympde/expr/evaluation.py
+++ b/sympde/expr/evaluation.py
@@ -31,8 +31,10 @@
 from sympde.topology.basic   import Boundary, Interface
 from sympde.topology.basic   import InteriorDomain
 from sympde.topology.domain  import NormalVector, TangentVector, NCube, NCubeInterior
-from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, MultiPatchMapping, JacobianInverseSymbol
-from sympde.topology.mapping import LogicalExpr, PullBack
+from sympde.topology.base_mapping import JacobianSymbol, JacobianInverseSymbol
+from sympde.topology.base_mapping import LogicalExpr, PullBack
+from sympde.topology.base_mapping import InterfaceMapping, MultiPatchMapping
+
 
 # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping
 
diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py
index efca3a98..5ef1856f 100644
--- a/sympde/expr/tests/test_expr.py
+++ b/sympde/expr/tests/test_expr.py
@@ -16,7 +16,7 @@
 
 from sympde.topology import dx1, dx2, dx3
 from sympde.topology import dx, dy, dz
-from sympde.topology import Mapping
+from sympde.topology import BaseMapping
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
 from sympde.topology import element_of, elements_of
 from sympde.topology import InteriorDomain, Union
@@ -867,7 +867,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      = BaseMapping('M', dim=3)
 
     mapped_domain = M(domain)
 
diff --git a/sympde/expr/tests/test_tensor_2d.py b/sympde/expr/tests/test_tensor_2d.py
index 0fc9e31d..f68e8667 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 BaseMapping
 
 from sympde.expr.expr import BilinearForm, integral
 from sympde.expr.evaluation import TensorExpr
@@ -59,11 +59,11 @@ def test_tensorize_2d_2():
     # ...
 
 #==============================================================================
-def test_tensorize_2d_1_mapping():
+def test_tensorize_2d_1_BaseMapping():
 
     DIM = 2
 
-    M = Mapping('Map', dim=DIM)
+    M = BaseMapping('Map', dim=DIM)
 
     domain = M(Domain('Omega', dim=DIM))
 
@@ -85,10 +85,10 @@ def test_tensorize_2d_1_mapping():
     # ...
 
 #==============================================================================
-def test_tensorize_2d_2_mapping():
+def test_tensorize_2d_2_BaseMapping():
 
     DIM = 2
-    M = Mapping('M', dim=DIM)
+    M = BaseMapping('M', dim=DIM)
     domain = M(Domain('Omega', dim=DIM))
 
     V = VectorFunctionSpace('V', domain)
diff --git a/sympde/topology/__init__.py b/sympde/topology/__init__.py
index de44bd3c..9b7eefce 100644
--- a/sympde/topology/__init__.py
+++ b/sympde/topology/__init__.py
@@ -1,9 +1,9 @@
-from .basic              import *
-from .derivatives        import *
-from .datatype           import *
-from .domain             import *
-from .mapping            import *
-from .measure            import *
-from .space              import *
-from .analytical_mapping import *
-from .callable_mapping   import *
+from .basic                     import *
+from .derivatives               import *
+from .datatype                  import *
+from .domain                    import *
+from .base_mapping              import *
+from .base_analytic_mapping     import *
+from .analytic_mappings         import *
+from .measure                   import *
+from .space                     import *
diff --git a/sympde/topology/analytical_mapping.py b/sympde/topology/analytic_mappings.py
similarity index 69%
rename from sympde/topology/analytical_mapping.py
rename to sympde/topology/analytic_mappings.py
index 05b9760b..ba998915 100644
--- a/sympde/topology/analytical_mapping.py
+++ b/sympde/topology/analytic_mappings.py
@@ -1,8 +1,8 @@
-from .mapping import Mapping
+from .base_analytic_mapping import BaseAnalyticMapping
 
-class IdentityMapping(Mapping):
+class IdentityMapping(BaseAnalyticMapping):
     """
-    Represents an identity 1D/2D/3D Mapping object.
+    Represents an identity 1D/2D/3D BaseAnalyticMapping object.
 
     Examples
 
@@ -12,9 +12,9 @@ class IdentityMapping(Mapping):
                     'z': 'x3'}
 
 #==============================================================================
-class AffineMapping(Mapping):
+class AffineMapping(BaseAnalyticMapping):
     """
-    Represents a 1D/2D/3D Affine Mapping object.
+    Represents a 1D/2D/3D Affine BaseAnalyticMapping object.
 
     Examples
 
@@ -24,9 +24,9 @@ class AffineMapping(Mapping):
                     'z': 'c3 + a31*x1 + a32*x2 + a33*x3'}
 
 #==============================================================================
-class PolarMapping(Mapping):
+class PolarMapping(BaseAnalyticMapping):
     """
-    Represents a Polar 2D Mapping object (Annulus).
+    Represents a Polar 2D BaseAnalyticMapping object (Annulus).
 
     Examples
 
@@ -38,9 +38,9 @@ class PolarMapping(Mapping):
     _pdim        = 2
 
 #==============================================================================
-class TargetMapping(Mapping):
+class TargetMapping(BaseAnalyticMapping):
     """
-    Represents a Target 2D Mapping object.
+    Represents a Target 2D BaseAnalyticMapping object.
 
     Examples
 
@@ -52,9 +52,9 @@ class TargetMapping(Mapping):
     _pdim        = 2
 
 #==============================================================================
-class CzarnyMapping(Mapping):
+class CzarnyMapping(BaseAnalyticMapping):
     """
-    Represents a Czarny 2D Mapping object.
+    Represents a Czarny 2D BaseAnalyticMapping object.
 
     Examples
 
@@ -67,19 +67,27 @@ class CzarnyMapping(Mapping):
     _pdim        = 2
 
 #==============================================================================
-class CollelaMapping2D(Mapping):
-    """
-    Represents a Collela 2D Mapping object.
+class CollelaMapping2D(BaseAnalyticMapping):
 
-    """
     _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.'}
 
-    _ldim        = 2
-    _pdim        = 2
+    _ldim = 2
+    _pdim = 2
+
+
+#==============================================================================
+class CollelaMapping3D(BaseAnalyticMapping):
+
+    _expressions = {'x': 'k1*(x1 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))',
+                    'y': 'k2*(x2 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))',
+                    'z': 'k3*x3'}
+
+    _ldim        = 3
+    _pdim        = 3
 
 #==============================================================================
-class TorusMapping(Mapping):
+class TorusMapping(BaseAnalyticMapping):
     """
     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 +106,7 @@ class TorusMapping(Mapping):
 
 #==============================================================================
 # TODO [YG, 07.10.2022]: add test in sympde/topology/tests/test_logical_expr.py
-class TorusSurfaceMapping(Mapping):
+class TorusSurfaceMapping(BaseAnalyticMapping):
     """
     3D surface obtained by "slicing" the torus above at r = a.
     The parametrization uses the coordinates (x1, x2) = (theta, phi), where:
@@ -116,7 +124,7 @@ class TorusSurfaceMapping(Mapping):
 
 #==============================================================================
 # TODO [YG, 07.10.2022]: add test in sympde/topology/tests/test_logical_expr.py
-class TwistedTargetSurfaceMapping(Mapping):
+class TwistedTargetSurfaceMapping(BaseAnalyticMapping):
     """
     3D surface obtained by "twisting" the TargetMapping out of the (x, y) plane
 
@@ -129,7 +137,7 @@ class TwistedTargetSurfaceMapping(Mapping):
     _pdim        = 3
 
 #==============================================================================
-class TwistedTargetMapping(Mapping):
+class TwistedTargetMapping(BaseAnalyticMapping):
     """
     3D volume obtained by "extruding" the TwistedTargetSurfaceMapping along z.
 
@@ -142,7 +150,7 @@ class TwistedTargetMapping(Mapping):
     _pdim        = 3
 
 #==============================================================================
-class SphericalMapping(Mapping):
+class SphericalMapping(BaseAnalyticMapping):
     """
     Parametrization of a sphere (or a portion of it) using spherical
     coordinates (x1, x2, x3) = (r, theta, phi), where:
@@ -158,3 +166,24 @@ class SphericalMapping(Mapping):
 
     _ldim        = 3
     _pdim        = 3
+
+#==============================================================================
+class Collela3D(BaseAnalyticMapping):
+
+    _expressions = {'x':'2.*(x1 + 0.1*sin(2.*pi*x1)*sin(2.*pi*x2)) - 1.',
+                    'y':'2.*(x2 + 0.1*sin(2.*pi*x1)*sin(2.*pi*x2)) - 1.',
+                    'z':'2.*x3  - 1.'}
+
+#==============================================================================
+class TransposedPolarMapping(BaseAnalyticMapping):
+    """
+    Represents a Transposed (x1 <> x2) Polar 2D Mapping object (Annulus).
+
+    Examples
+
+    """
+    _expressions = {'x': 'c1 + (rmin*(1-x2)+rmax*x2)*cos(x1)',
+                    'y': 'c2 + (rmin*(1-x2)+rmax*x2)*sin(x1)'}
+
+    _ldim        = 2
+    _pdim        = 2
diff --git a/sympde/topology/base_analytic_mapping.py b/sympde/topology/base_analytic_mapping.py
new file mode 100644
index 00000000..55c49710
--- /dev/null
+++ b/sympde/topology/base_analytic_mapping.py
@@ -0,0 +1,216 @@
+# coding: utf-8
+
+import numpy as np
+import itertools as it
+
+from sympy                 import lambdify
+from sympy.core            import Symbol
+
+from .basic       import BasicDomain
+from .base_mapping import BaseMapping, MappedDomain
+
+
+# TODO fix circular dependency between sympde.topology.domain and sympde.topology.mapping
+# TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping
+
+__all__ = (
+    'lambdify_sympde',
+    'BaseAnalyticMapping',
+)
+
+
+def lambdify_sympde(variables, expr):
+    """
+    Custom lambify function that covers the
+    shortcomings of sympy's lambdify. Most notably,
+    this function uses numpy broadcasting rules to
+    compute the shape of the output.
+
+    Parameters
+    ----------
+    variables : sympy.core.symbol.Symbol or list of sympy.core.symbol.Symbol
+        variables that appear in the expression
+    expr :
+        Sympy expression
+
+    Returns
+    -------
+    lambda_f : callable
+        Lambdified function built using numpy.
+
+    Notes
+    -----
+    Compared to Sympy's lambdify, this function
+    is capable of properly handling constant values,
+    and array_like structures where not all components
+    depend on all variables. See below.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from sympy import symbols,  Matrix
+    >>> from sympde.utilities.utils import lambdify_sympde
+    >>> x, y = symbols("x,y")
+    >>> expr = Matrix([[x, x + y], [0, y]])
+    >>> f = lambdify_sympde([x,y], expr)
+    >>> f(np.array([[0, 1]]), np.array([[2], [3]]))
+    array([[[[0., 1.],
+             [0., 1.]],
+
+            [[2., 3.],
+             [3., 4.]]],
+
+
+           [[[0., 0.],
+             [0., 0.]],
+
+            [[2., 2.],
+             [3., 3.]]]])
+    """
+    array_expr = np.asarray(expr)
+    scalar_shape = array_expr.shape
+    if scalar_shape == ():
+        f = lambdify(variables, expr, 'numpy')
+        def f_vec_sc(*XYZ):
+            b = np.broadcast(*XYZ)
+            if b.ndim == 0:
+                return f(*XYZ)
+            temp = np.asarray(f(*XYZ))
+            if b.shape == temp.shape:
+                return temp
+
+            result = np.zeros(b.shape)
+            result[...] = temp
+            return result
+        return f_vec_sc
+
+    else:
+        scalar_functions = {}
+        for multi_index in it.product(*tuple(range(s) for s in scalar_shape)):
+            scalar_functions[multi_index] = lambdify(variables, array_expr[multi_index], 'numpy')
+
+        def f_vec_v(*XYZ):
+            b = np.broadcast(*XYZ)
+            result = np.zeros(scalar_shape + b.shape)
+            for multi_index in it.product(*tuple(range(s) for s in scalar_shape)):
+                result[multi_index] = scalar_functions[multi_index](*XYZ)
+            return result
+        return f_vec_v
+
+
+#==============================================================================
+class BaseAnalyticMapping(BaseMapping):
+    """
+    Represents a BaseAnalyticMapping object.
+
+    Examples
+
+    """
+
+    def __new__(cls, name, dim=None, **kwargs):
+
+        obj = super().__new__(cls, name, dim, **kwargs)
+
+        if obj.expressions :
+            obj._func_eval = tuple(lambdify_sympde( obj._logical_coordinates, expr) for expr in obj._expressions)
+            obj._jac_eval = lambdify_sympde( obj._logical_coordinates, obj._jac)
+            obj._inv_jac_eval = lambdify_sympde( obj._logical_coordinates, obj._inv_jac)
+            obj._metric_eval = lambdify_sympde( obj._logical_coordinates, obj._metric)
+            obj._metric_det_eval = lambdify_sympde( obj._logical_coordinates, obj._metric_det)
+        else:
+            raise TypeError("BaseAnalyticMapping should have an expression")
+
+        return obj
+
+    #--------------------------------------------------------------------------
+    #Abstract Interface :
+
+    def _evaluate_domain( self, domain ):
+        assert(isinstance(domain, BasicDomain))
+        return MappedDomain(self, domain)
+
+    def _evaluate( self, *Xs ):
+        #int, float or numpy arrays
+        if self._func_eval is None :
+            raise TypeError("not a callable object")
+        else :
+            assert len(Xs)==self.ldim
+            Xshape = np.shape(Xs[0])
+            for X in Xs:
+                assert np.shape(X) == Xshape
+            return tuple( f( *Xs ) for f in self._func_eval)
+
+    def _jacobian_evaluate( self, *Xs ):
+        #int, float or numpy arrays
+        if self._jac_eval is None:
+            raise TypeError("not a callable object")
+        else :
+            assert len(Xs)==self.ldim
+            Xshape = np.shape(Xs[0])
+            for X in Xs:
+                assert np.shape(X) == Xshape
+            return self._jac_eval(*Xs)
+
+    def _jacobian_inv_evaluate( self, *Xs ):
+        #int, float or numpy arrays
+        if self._inv_jac_eval is None:
+            raise TypeError("not a callable object")
+        else :
+            assert len(Xs)==self.ldim
+            Xshape = np.shape(Xs[0])
+            for X in Xs:
+                assert np.shape(X) == Xshape
+            return  self._inv_jac_eval(*Xs)
+
+    def _metric_evaluate( self, *Xs ):
+        if self._metric_eval is None:
+            raise TypeError("not a callable object")
+        else :
+            assert len(Xs)==self.ldim
+            Xshape = np.shape(Xs[0])
+            for X in Xs:
+                assert np.shape(X) == Xshape
+            return  self._metric_eval(*Xs)
+
+    def _metric_det_evaluate( self, *Xs ):
+        if self._metric_det_eval is None:
+            raise TypeError("not a callable object")
+        else :
+            assert len(Xs)==self.ldim
+            Xshape = np.shape(Xs[0])
+            for X in Xs:
+                assert np.shape(X) == Xshape
+            return self._metric_det_eval(*Xs)
+
+    def __call__( self, *args ):
+        if len(args) == 1 and isinstance(args[0], BasicDomain):
+            return self._evaluate_domain(args[0])
+        elif all(isinstance(arg, (int, float, Symbol, np.ndarray)) for arg in args):
+            return self._evaluate(*args)
+        else:
+            raise TypeError("Invalid arguments for __call__")
+
+    def jacobian_eval( self, *args ):
+        if all(isinstance(arg, (int, float, Symbol, np.ndarray)) for arg in args):
+            return self._jacobian_evaluate(*args)
+        else:
+            raise TypeError("Invalid arguments for jacobian_eval")
+
+    def jacobian_inv_eval( self, *args ):
+        if all(isinstance(arg, (int, float, Symbol, np.ndarray)) for arg in args):
+            return self._jacobian_inv_evaluate(*args)
+        else:
+            raise TypeError("Invalid arguments for jacobian_inv_eval")
+
+    def metric_eval( self, *args ):
+        if all(isinstance(arg, (int, float, Symbol, np.ndarray)) for arg in args):
+            return self._metric_evaluate(*args)
+        else:
+            raise TypeError("Invalid arguments for metric_eval")
+
+    def metric_det_eval( self, *args ):
+        if all(isinstance(arg, (int, float, Symbol, np.ndarray)) for arg in args):
+            return self._metric_det_evaluate(*args)
+        else:
+            raise TypeError("Invalid arguments for metric_det_eval")
+
diff --git a/sympde/topology/mapping.py b/sympde/topology/base_mapping.py
similarity index 92%
rename from sympde/topology/mapping.py
rename to sympde/topology/base_mapping.py
index 76090fc0..c77a788d 100644
--- a/sympde/topology/mapping.py
+++ b/sympde/topology/base_mapping.py
@@ -14,9 +14,9 @@
 from sympy                 import sqrt, symbols
 from sympy.core.exprtools  import factor_terms
 from sympy.polys.polytools import parallel_poly_from_expr
+import numpy as np
 
 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.calculus.core     import PlusInterfaceOperator, MinusInterfaceOperator
@@ -40,11 +40,12 @@
 from .derivatives import get_atom_logical_derivatives, get_index_logical_derivatives_atom
 from .derivatives import LogicalGrad_1d, LogicalGrad_2d, LogicalGrad_3d
 
+
 # TODO fix circular dependency between sympde.topology.domain and sympde.topology.mapping
 # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping
 
 __all__ = (
-    'BasicCallableMapping',
+    'BaseMapping',
     'Contravariant',
     'Covariant',
     'InterfaceMapping',
@@ -60,6 +61,7 @@
     'SymbolicExpr',
     'SymbolicWeightedVolume',
     'get_logical_test_function',
+    'numpy_to_native_python',
 )
 
 #==============================================================================
@@ -84,63 +86,15 @@ def get_logical_test_function(u):
     el              = l_space.element(u.name)
     return el
 
-import numpy as np
-
 def numpy_to_native_python(a):
     if isinstance(a, np.generic):
         return a.item()
     return a
 
 #==============================================================================
-class BasicCallableMapping(ABC):
+class BaseMapping(IndexedBase):
     """
-    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):
-    """
-    Represents a Mapping object.
+    Represents a BaseMapping object.
 
     Examples
 
@@ -149,7 +103,6 @@ class Mapping(BasicMapping):
     _jac          = None
     _inv_jac      = None
     _constants    = None
-    _callable_map = None
     _ldim         = None
     _pdim         = None
 
@@ -177,7 +130,6 @@ def __new__(cls, name, dim=None, **kwargs):
             ldim = dim
             pdim = dim
 
-
         obj = IndexedBase.__new__(cls, name, shape=pdim)
 
         if not evaluate:
@@ -262,29 +214,6 @@ 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
@@ -372,7 +301,7 @@ def set_plus_minus( self, **kwargs):
         self._is_minus = minus
 
     def copy(self):
-        obj = Mapping(self.name,
+        obj = BaseMapping(self.name,
                      ldim=self.ldim,
                      pdim=self.pdim,
                      evaluate=False)
@@ -389,7 +318,6 @@ def copy(self):
         obj._inv_jac             = self._inv_jac
         obj._metric              = self._metric
         obj._metric_det          = self._metric_det
-        obj.__callable_map       = self._callable_map
         obj._is_plus             = self._is_plus
         obj._is_minus            = self._is_minus
         return obj
@@ -407,21 +335,21 @@ def _sympystr(self, printer):
         return sstr(self.name)
 
 #==============================================================================
-class InverseMapping(Mapping):
+class InverseMapping(BaseMapping):
     def __new__(cls, mapping):
-        assert isinstance(mapping, Mapping)
+        assert isinstance(mapping, BaseMapping)
         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)
+        return BaseMapping.__new__(cls, name, ldim=ldim, pdim=pdim, coordinates=coords, jacobian=jacobian)
 
 #==============================================================================
 class JacobianSymbol(MatrixSymbolicExpr):
     _axis = None
     def __new__(cls, mapping, axis=None):
-        assert isinstance(mapping, Mapping)
+        assert isinstance(mapping, BaseMapping)
         if axis is not None:
             assert isinstance(axis, (int, Integer))
         obj = MatrixSymbolicExpr.__new__(cls, mapping)
@@ -449,7 +377,7 @@ def __hash__(self):
         return hash(self._hashable_content())
 
     def _eval_subs(self, old, new):
-        if isinstance(new, Mapping):
+        if isinstance(new, BaseMapping):
             if self.axis is not None:
                 obj = JacobianSymbol(new, self.axis)
             else:
@@ -468,7 +396,7 @@ class JacobianInverseSymbol(MatrixSymbolicExpr):
     _axis = None
     is_Matrix     = False
     def __new__(cls, mapping, axis=None):
-        assert isinstance(mapping, Mapping)
+        assert isinstance(mapping, BaseMapping)
         if axis is not None:
             assert isinstance(axis, int)
         obj = MatrixSymbolicExpr.__new__(cls, mapping)
@@ -500,21 +428,21 @@ def _sympystr(self, printer):
             return 'Jacobian({})**(-1)'.format(sstr(self.mapping.name))
 
 #==============================================================================
-class InterfaceMapping(Mapping):
+class InterfaceMapping(BaseMapping):
     """
     InterfaceMapping is used to represent a mapping in the interface.
 
     Attributes
     ----------
-    minus : Mapping
+    minus : BaseMapping
         the mapping on the negative direction of the interface
-    plus  : Mapping
+    plus  : BaseMapping
         the mapping on the positive direction of the interface
     """
 
     def __new__(cls, minus, plus):
-        assert isinstance(minus, Mapping)
-        assert isinstance(plus,  Mapping)
+        assert isinstance(minus, BaseMapping)
+        assert isinstance(plus,  BaseMapping)
         minus = minus.copy()
         plus  = plus.copy()
 
@@ -522,7 +450,7 @@ def __new__(cls, minus, plus):
         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  = BaseMapping.__new__(cls, name, ldim=minus.ldim, pdim=minus.pdim)
         obj._minus = minus
         obj._plus  = plus
         return obj
@@ -548,7 +476,7 @@ def _eval_simplify(self, **kwargs):
         return self
 
 #==============================================================================
-class MultiPatchMapping(Mapping):
+class MultiPatchMapping(BaseMapping):
 
     def __new__(cls, dic):
         assert isinstance( dic, dict)
@@ -570,10 +498,6 @@ def ldim(self):
     def pdim(self):
         return list(self.mappings.values())[0].pdim
 
-    @property
-    def is_analytical(self):
-        return all(e.is_analytical for e in self.mappings.values())
-
     def _eval_subs(self, old, new):
         return self
 
@@ -594,7 +518,7 @@ class MappedDomain(BasicDomain):
 
     @cacheit
     def __new__(cls, mapping, logical_domain):
-        assert(isinstance(mapping, Mapping))
+        assert(isinstance(mapping, BaseMapping))
         assert(isinstance(logical_domain, BasicDomain))
         if isinstance(logical_domain, Domain):
             kwargs = dict(
@@ -744,7 +668,7 @@ def eval(cls, F):
 
         Parameters:
         ----------
-         F: Mapping
+         F: BaseMapping
             mapping object
 
         Returns:
@@ -753,8 +677,8 @@ def eval(cls, F):
             the jacobian matrix
         """
 
-        if not isinstance(F, Mapping):
-            raise TypeError('> Expecting a Mapping object')
+        if not isinstance(F, BaseMapping):
+            raise TypeError('> Expecting a BaseMapping object')
 
         if F.jacobian_expr is not None:
             return F.jacobian_expr
@@ -792,7 +716,7 @@ def eval(cls, F, v):
 
         Parameters:
         ----------
-         F: Mapping
+         F: BaseMapping
             mapping object
 
          v: <tuple|list|Tuple|ImmutableDenseMatrix|Matrix>
@@ -841,7 +765,7 @@ def eval(cls, F, v):
 
         Parameters:
         ----------
-         F: Mapping
+         F: BaseMapping
             mapping object
 
          v: <tuple|list|Tuple|ImmutableDenseMatrix|Matrix>
@@ -853,8 +777,8 @@ def eval(cls, F, v):
             the contravariant transformation
         """
 
-        if not isinstance(F, Mapping):
-            raise TypeError('> Expecting a Mapping')
+        if not isinstance(F, BaseMapping):
+            raise TypeError('> Expecting a BaseMapping')
 
         if not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)):
             raise TypeError('> Expecting a tuple, list, Tuple, Matrix')
@@ -1207,7 +1131,6 @@ def eval(cls, expr, domain, **options):
             domain  = expr.domain
             mapping = domain.mapping
 
-
             assert domain is not None
 
             if expr.is_domain_integral:
@@ -1346,7 +1269,7 @@ def eval(cls, *_args, **kwargs):
 
         elif isinstance(expr, Indexed):
             base = expr.base
-            if isinstance(base, Mapping):
+            if isinstance(base, BaseMapping):
                 if expr.indices[0] == 0:
                     name = 'x'
                 elif expr.indices[0] == 1:
@@ -1391,7 +1314,7 @@ def eval(cls, *_args, **kwargs):
                     code += k*n
             return cls.eval(atom, code=code)
 
-        elif isinstance(expr, Mapping):
+        elif isinstance(expr, BaseMapping):
             return Symbol(expr.name)
 
         # ... this must be done here, otherwise codegen for FEM will not work
diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py
index 07050c61..e770defc 100644
--- a/sympde/topology/basic.py
+++ b/sympde/topology/basic.py
@@ -324,7 +324,7 @@ def rotate(self, *directions):
         raise NotImplementedError('only 2d case is available')
 
     def join(self, boundary, ornt=None):
-        from sympde.topology.mapping import InterfaceMapping
+        from sympde.topology.base_mapping import InterfaceMapping
         # TODO be careful with '|' in psydac
         if self.mapping and boundary.mapping:
             int_map            = InterfaceMapping(self.mapping , boundary.mapping)
diff --git a/sympde/topology/callable_mapping.py b/sympde/topology/callable_mapping.py
deleted file mode 100644
index f1f4d60f..00000000
--- a/sympde/topology/callable_mapping.py
+++ /dev/null
@@ -1,95 +0,0 @@
-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
-
diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py
index f5d3537a..4a1f5c42 100644
--- a/sympde/topology/derivatives.py
+++ b/sympde/topology/derivatives.py
@@ -29,6 +29,7 @@
 
 from .space   import ScalarFunction, VectorFunction, IndexedVectorFunction
 
+
 #==============================================================================
 class DifferentialOperator(LinearOperator):
     """
@@ -51,7 +52,6 @@ def __getitem__(self, indices, **kw_args):
     @classmethod
     @cacheit
     def eval(cls, expr):
-
         types = (VectorFunction, ScalarFunction, DifferentialOperator)
 
         if isinstance(expr, _logical_partial_derivatives):
@@ -97,7 +97,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, BaseMapping):
             return cls(expr, evaluate=False)
         elif not has(expr, types):
             if expr.is_number:
@@ -106,7 +106,7 @@ def eval(cls, expr):
             elif isinstance(expr, Expr):
                 x = Symbol(cls.coordinate)
                 if cls.logical:
-                    M = expr.atoms(Mapping)
+                    M = expr.atoms(BaseMapping)
                     if len(M)>0:
                         M = list(M)[0]
                         expr_primes = [diff(expr, M[i]) for i in range(M.pdim)]
@@ -1288,5 +1288,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 .base_mapping import BaseMapping, Jacobian
diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py
index 1acf9a5d..79b104f5 100644
--- a/sympde/topology/domain.py
+++ b/sympde/topology/domain.py
@@ -23,10 +23,9 @@
 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.base_mapping import BaseMapping
 # TODO add pdim
 
 iterable_types = (tuple, list, Tuple, Union)
@@ -46,7 +45,7 @@ def __new__(cls, name : str, *,
             boundaries : TypeUnion[Iterable[Boundary], Boundary, None] = None,
             dim : Optional[int] = None,
             connectivity : Optional[Connectivity] = None,
-            mapping : Optional[Mapping] = None,
+            mapping : Optional[BaseMapping] = None,
             logical_domain : Optional[Domain] = None):
         """
         Interiors or connectivity must be given. When the mapping is given 
@@ -157,6 +156,7 @@ def __new__(cls, name : str, *,
         obj._dtype          = dtype
         obj._dim            = dim
         obj._logical_domain = logical_domain
+        obj._mapping = mapping
         return obj
 
     @property
@@ -176,9 +176,14 @@ def boundary(self) -> TypeUnion[Union, Boundary]:
         return self.args[2]
 
     @property
-    def mapping(self) -> Optional[Mapping]:
+    def mapping(self) -> Optional[BaseMapping]:
         """The mapping that maps the logical domain to the physical domain"""
-        return self.args[3]
+        return self._mapping
+
+    @mapping.setter
+    def mapping(self, value: BaseMapping):
+        """Set the mapping that maps the logical domain to the physical domain"""
+        self._mapping = value
 
     @property
     def logical_domain(self) -> Domain:
@@ -337,7 +342,7 @@ def from_file( cls, filename ):
         if not(ext == '.h5'):
             raise ValueError('> Only h5 files are supported')
         # ...
-        from sympde.topology.mapping import Mapping, MultiPatchMapping
+        from sympde.topology.base_mapping import BaseMapping
 
         h5  = h5py.File( filename, mode='r' )
         yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader )
@@ -359,7 +364,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     = [BaseMapping(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)}
 
@@ -486,7 +491,6 @@ def join(cls, patches, connectivity, name):
 
         patch_given_by_indices = (len(connectivity) > 0 and isinstance(connectivity[0][0][0], int))
 
-        from sympde.topology.mapping import MultiPatchMapping
         # ... connectivity
         interfaces = {}
         boundaries = []
@@ -536,6 +540,7 @@ def join(cls, patches, connectivity, name):
             for k,v in connectivity.items():
                 logical_connectivity[v.logical_domain.name] = v.logical_domain
 
+            from sympde.topology.base_mapping import MultiPatchMapping
             mapping        = MultiPatchMapping({e.logical_domain: e.mapping for e in interiors})
             logical_domain = Domain(name,
                             interiors=logical_interiors,
diff --git a/sympde/topology/tests/test_mapping.py b/sympde/topology/tests/test_base_mapping.py
similarity index 92%
rename from sympde/topology/tests/test_mapping.py
rename to sympde/topology/tests/test_base_mapping.py
index 8d7c6798..e192d32d 100644
--- a/sympde/topology/tests/test_mapping.py
+++ b/sympde/topology/tests/test_base_mapping.py
@@ -7,19 +7,20 @@
 from sympy.tensor import IndexedBase
 from sympy import symbols, simplify
 
-from sympde.topology import Mapping, MappedDomain, AffineMapping
-from sympde.topology import dx, dy, dz
-from sympde.topology import dx1, dx2, dx3
-from sympde.topology import Domain
+from sympde.topology.domain            import Domain
+from sympde.topology.base_mapping      import BaseMapping, MappedDomain
+from sympde.topology.analytic_mappings import AffineMapping
+from sympde.topology.derivatives       import dx, dy, dz
+from sympde.topology.derivatives       import dx1, dx2, dx3
 
-from sympde.topology.mapping import Jacobian, Covariant, Contravariant
+from sympde.topology.base_mapping import Jacobian, Covariant, Contravariant
 # ...
 def test_mapping_1d():
     print('============ test_mapping_1d ==============')
 
     dim = 1
 
-    F = Mapping('F', dim=dim)
+    F = BaseMapping('F', dim=dim)
 
     assert(F.name == 'F')
 
@@ -40,7 +41,7 @@ def test_mapping_2d():
 
     dim = 2
 
-    F = Mapping('F', dim=dim)
+    F = BaseMapping('F', dim=dim)
 
     a,b = symbols('a b')
     ab = Tuple(a, b)
@@ -80,7 +81,7 @@ def test_mapping_3d():
 
     dim = 3
 
-    F = Mapping('F', dim=dim)
+    F = BaseMapping('F', dim=dim)
 
     a,b,c = symbols('a b c')
     abc = Tuple(a, b, c)
@@ -125,10 +126,11 @@ def test_mapping_2d_2():
     print('============ test_mapping_2d_2 ==============')
 
     dim   = 2
-    F      = Mapping('F', dim=dim)
+    F      = BaseMapping('F', dim=dim)
     domain = Domain('Omega', dim=dim)
     D      = F(domain)
 
+# ...
 def test_AffineMapping():
     print('============ test_AffineMapping ==============')
 
@@ -158,3 +160,9 @@ def teardown_function():
     from sympy.core import cache
     cache.clear_cache()
 
+if __name__ == '__main__' :
+    test_mapping_1d()
+    test_mapping_2d()
+    test_mapping_3d()
+    test_mapping_2d_2()
+    test_AffineMapping()
diff --git a/sympde/topology/tests/test_derivatives.py b/sympde/topology/tests/test_derivatives.py
index 3694adfe..3f88d112 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 BaseMapping
 
 
 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      = BaseMapping('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      = BaseMapping('M', dim=2)
 
     mapped_domain = M(domain)
 
diff --git a/sympde/topology/tests/test_logical_expr.py b/sympde/topology/tests/test_logical_expr.py
index 99108c4d..7825676d 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, BaseMapping, Square
 from sympde.topology import dx, dy
 from sympde.topology import dx1, dx2, dx3
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
@@ -27,13 +27,13 @@
 from sympde.expr     import BilinearForm, integral
 from sympde.calculus import grad, div, curl, dot
 
-from sympde.topology.mapping import Jacobian
+from sympde.topology.base_mapping import Jacobian
 
 #==============================================================================
 def test_logical_expr_1d_1():
     dim = 1
 
-    M = Mapping('M', dim=dim)
+    M = BaseMapping('M', dim=dim)
     domain = M(Domain('Omega', dim=dim))
 
     alpha = Constant('alpha')
@@ -80,7 +80,7 @@ def test_logical_expr_1d_1():
 def test_symbolic_expr_1d_1():
     dim = 1
 
-    M = Mapping('M', dim=dim)
+    M = BaseMapping('M', dim=dim)
     domain = M(Domain('Omega', dim=dim))
 
     V = ScalarFunctionSpace('V', domain, kind='h1')
@@ -149,7 +149,7 @@ def test_symbolic_expr_1d_1():
 def test_logical_expr_2d_1():
     dim = 2
 
-    M = Mapping('M', dim=dim)
+    M = BaseMapping('M', dim=dim)
     domain = M(Domain('Omega', dim=dim))
 
     alpha = Constant('alpha')
@@ -215,8 +215,8 @@ def test_logical_expr_2d_2():
     A = Square('A')
     B = Square('B')
 
-    M1 = Mapping('M1', dim=dim)
-    M2 = Mapping('M2', dim=dim)
+    M1 = BaseMapping('M1', dim=dim)
+    M2 = BaseMapping('M2', dim=dim)
 
     D1 = M1(A)
     D2 = M2(B)
@@ -268,8 +268,8 @@ def test_logical_expr_2d_3():
     A = Square('A')
     B = Square('B')
 
-    M1 = Mapping('M1', dim=dim)
-    M2 = Mapping('M2', dim=dim)
+    M1 = BaseMapping('M1', dim=dim)
+    M2 = BaseMapping('M2', dim=dim)
 
     D1 = M1(A)
     D2 = M2(B)
@@ -292,7 +292,7 @@ def test_logical_expr_2d_3():
 #==============================================================================
 def test_symbolic_expr_2d_1():
     dim    = 2
-    M      = Mapping('M', dim=dim)
+    M      = BaseMapping('M', dim=dim)
     domain = M(Domain('Omega', dim=dim))
 
     V = ScalarFunctionSpace('V', domain, kind='h1')
@@ -368,7 +368,7 @@ def test_symbolic_expr_2d_1():
 def test_logical_expr_3d_1():
     dim = 3
 
-    M = Mapping('M', dim=dim)
+    M = BaseMapping('M', dim=dim)
     domain = M(Domain('Omega', dim=dim))
 
     alpha = Constant('alpha')
@@ -421,7 +421,7 @@ def test_logical_expr_3d_2():
 
     dim   = 3
     domain = Domain('Omega', dim=dim)
-    M      = Mapping('M', dim=dim)
+    M      = BaseMapping('M', dim=dim)
 
     mapped_domain = M(domain)
 
@@ -443,7 +443,7 @@ def test_logical_expr_3d_3():
 
     dim    = 3
     domain = Domain('Omega', dim=dim)
-    M      = Mapping('M', dim=dim)
+    M      = BaseMapping('M', dim=dim)
 
     mapped_domain = M(domain)
 
@@ -464,7 +464,7 @@ def test_logical_expr_3d_4():
 
     dim    = 3
     domain = Domain('Omega', dim=dim)
-    M      = Mapping('M', dim=dim)
+    M      = BaseMapping('M', dim=dim)
 
     mapped_domain = M(domain)
 
@@ -485,7 +485,7 @@ def test_logical_expr_3d_5():
 
     dim    = 3
     domain = Domain('Omega', dim=dim)
-    M      = Mapping('M', dim=dim)
+    M      = BaseMapping('M', dim=dim)
 
     mapped_domain = M(domain)
 
@@ -509,7 +509,7 @@ def test_logical_expr_3d_5():
 #==============================================================================
 def test_symbolic_expr_3d_1():
     dim    = 3
-    M      = Mapping('M', dim=dim)
+    M      = BaseMapping('M', dim=dim)
     domain = M(Domain('Omega', dim=dim))
 
     V = ScalarFunctionSpace('V', domain, kind='h1')
diff --git a/sympde/topology/tests/test_callable_mapping.py b/sympde/topology/tests/test_mapping_evaluation.py
similarity index 60%
rename from sympde/topology/tests/test_callable_mapping.py
rename to sympde/topology/tests/test_mapping_evaluation.py
index 23462c3d..535d342a 100644
--- a/sympde/topology/tests/test_callable_mapping.py
+++ b/sympde/topology/tests/test_mapping_evaluation.py
@@ -1,8 +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
+from sympde.topology.analytic_mappings import IdentityMapping, AffineMapping, PolarMapping
 
 # Tolerance for testing float equality
 RTOL = 1e-15
@@ -11,8 +9,7 @@
 #==============================================================================
 def test_identity_mapping_1d():
 
-    F = IdentityMapping('F', dim=1)
-    f = F.get_callable_mapping()
+    f = IdentityMapping('F', dim=1)
 
     assert f.ldim == 1
     assert f.pdim == 1
@@ -21,16 +18,15 @@ def test_identity_mapping_1d():
 
     for x1 in x1_pts:
         assert f(x1)              == (x1,)
-        assert f.jacobian(x1)     == 1.0
-        assert f.jacobian_inv(x1) == 1.0
-        assert f.metric(x1)       == 1.0
-        assert f.metric_det(x1)   == 1.0
+        assert f.jacobian_eval(x1)     == 1.0
+        assert f.jacobian_inv_eval(x1) == 1.0
+        assert f.metric_eval(x1)       == 1.0
+        assert f.metric_det_eval(x1)   == 1.0
 
 
 def test_identity_mapping_2d():
 
-    F = IdentityMapping('F', dim=2)
-    f = F.get_callable_mapping()
+    f = IdentityMapping('F', dim=2)
 
     assert f.ldim == 2
     assert f.pdim == 2
@@ -44,16 +40,15 @@ def test_identity_mapping_2d():
     for x1 in x1_pts:
         for x2 in x2_pts:
             assert f(x1, x2) == (x1, x2)
-            assert np.array_equal(f.jacobian    (x1, x2), I)
-            assert np.array_equal(f.jacobian_inv(x1, x2), I)
-            assert np.array_equal(f.metric      (x1, x2), I)
-            assert f.metric_det(x1, x2) == 1.0
+            assert np.array_equal(f.jacobian_eval    (x1, x2), I)
+            assert np.array_equal(f.jacobian_inv_eval(x1, x2), I)
+            assert np.array_equal(f.metric_eval      (x1, x2), I)
+            assert f.metric_det_eval(x1, x2) == 1.0
 
 
 def test_identity_mapping_3d():
 
-    F = IdentityMapping('F', dim=3)
-    f = F.get_callable_mapping()
+    f = IdentityMapping('F', dim=3)
 
     assert f.ldim == 3
     assert f.pdim == 3
@@ -70,10 +65,10 @@ def test_identity_mapping_3d():
         for x2 in x2_pts:
             for x3 in x3_pts:
                 assert f(x1, x2, x3) == (x1, x2, x3)
-                assert np.array_equal(f.jacobian    (x1, x2, x3), I)
-                assert np.array_equal(f.jacobian_inv(x1, x2, x3), I)
-                assert np.array_equal(f.metric      (x1, x2, x3), I)
-                assert f.metric_det(x1, x2, x3) == 1.0
+                assert np.array_equal(f.jacobian_eval    (x1, x2, x3), I)
+                assert np.array_equal(f.jacobian_inv_eval(x1, x2, x3), I)
+                assert np.array_equal(f.metric_eval      (x1, x2, x3), I)
+                assert f.metric_det_eval(x1, x2, x3) == 1.0
 
 
 #------------------------------------------------------------------------------
@@ -82,8 +77,7 @@ def test_affine_mapping_1d():
     # x = 1 - x1
     params = {'c1': 1, 'a11': -1}
 
-    F = AffineMapping('F', **params, dim=1)
-    f = F.get_callable_mapping()
+    f = AffineMapping('F', **params, dim=1)
 
     assert f.ldim == 1
     assert f.pdim == 1
@@ -93,10 +87,10 @@ def test_affine_mapping_1d():
     assert f(1  ) == (0  ,)
 
     for x1 in [0, 0.5, 1]:
-        assert f.jacobian(x1)     == -1
-        assert f.jacobian_inv(x1) == -1
-        assert f.metric(x1)       ==  1
-        assert f.metric_det(x1)   ==  1
+        assert f.jacobian_eval(x1)     == -1
+        assert f.jacobian_inv_eval(x1) == -1
+        assert f.metric_eval(x1)       ==  1
+        assert f.metric_det_eval(x1)   ==  1
 
 
 def test_affine_mapping_2d():
@@ -106,8 +100,7 @@ def test_affine_mapping_2d():
 
     params = dict(c1=c1, c2=c2, a11=a11, a12=a12, a21=a21, a22=a22)
 
-    F = AffineMapping('F', **params, dim=2)
-    f = F.get_callable_mapping()
+    f = AffineMapping('F', **params, dim=2)
 
     assert f.ldim == 2
     assert f.pdim == 2
@@ -128,10 +121,10 @@ def test_affine_mapping_2d():
             x = c1 + a11 * x1 + a12 * x2
             y = c2 + a21 * x1 + a22 * x2
             assert f(x1, x2) == (x, y)
-            assert np.array_equal(f.jacobian    (x1, x2), J    )
-            assert np.array_equal(f.jacobian_inv(x1, x2), J_inv)
-            assert np.array_equal(f.metric      (x1, x2), G    )
-            assert f.metric_det(x1, x2) == G_det
+            assert np.array_equal(f.jacobian_eval    (x1, x2), J    )
+            assert np.array_equal(f.jacobian_inv_eval(x1, x2), J_inv)
+            assert np.array_equal(f.metric_eval      (x1, x2), G    )
+            assert f.metric_det_eval(x1, x2) == G_det
 
 
 def test_affine_mapping_3d():
@@ -147,8 +140,7 @@ def test_affine_mapping_3d():
             a21=a21, a22=a22, a23=a23,
             a31=a31, a32=a32, a33=a33,)
 
-    F = AffineMapping('F', **params, dim=3)
-    f = F.get_callable_mapping()
+    f = AffineMapping('F', **params, dim=3)
 
     assert f.ldim == 3
     assert f.pdim == 3
@@ -178,10 +170,10 @@ def test_affine_mapping_3d():
                 y = c2 + a21 * x1 + a22 * x2 + a23 * x3
                 z = c3 + a31 * x1 + a32 * x2 + a33 * x3
                 assert np.allclose(f(x1, x2, x3), (x, y, z), rtol=RTOL, atol=ATOL)
-                assert np.array_equal(f.jacobian    (x1, x2, x3), J    )
-                assert np.array_equal(f.jacobian_inv(x1, x2, x3), J_inv)
-                assert np.array_equal(f.metric      (x1, x2, x3), G    )
-                assert f.metric_det(x1, x2, x3) == G_det
+                assert np.array_equal(f.jacobian_eval    (x1, x2, x3), J    )
+                assert np.array_equal(f.jacobian_inv_eval(x1, x2, x3), J_inv)
+                assert np.array_equal(f.metric_eval      (x1, x2, x3), G    )
+                assert f.metric_det_eval(x1, x2, x3) == G_det
 
 
 #------------------------------------------------------------------------------
@@ -192,8 +184,7 @@ def test_polar_mapping():
 
     params = dict(c1=c1, c2=c2, rmin=rmin, rmax=rmax)
 
-    F = PolarMapping('F', **params)
-    f = F.get_callable_mapping()
+    f = PolarMapping('F', **params)
 
     assert f.ldim == 2
     assert f.pdim == 2
@@ -221,10 +212,10 @@ def test_polar_mapping():
                      [-J[1][0] / J_det,  J[0][0] / J_det]]
 
             assert f(x1, x2) == (x, y)
-            assert np.allclose(f.jacobian    (x1, x2), J    , rtol=RTOL, atol=ATOL)
-            assert np.allclose(f.jacobian_inv(x1, x2), J_inv, rtol=RTOL, atol=ATOL)
-            assert np.allclose(f.metric      (x1, x2), G    , rtol=RTOL, atol=ATOL)
-            assert np.allclose(f.metric_det  (x1, x2), G_det, rtol=RTOL, atol=ATOL)
+            assert np.allclose(f.jacobian_eval    (x1, x2), J    , rtol=RTOL, atol=ATOL)
+            assert np.allclose(f.jacobian_inv_eval(x1, x2), J_inv, rtol=RTOL, atol=ATOL)
+            assert np.allclose(f.metric_eval      (x1, x2), G    , rtol=RTOL, atol=ATOL)
+            assert np.allclose(f.metric_det_eval  (x1, x2), G_det, rtol=RTOL, atol=ATOL)
 
 
 #------------------------------------------------------------------------------
@@ -232,8 +223,7 @@ def test_polar_mapping():
 
 def test_identity_mapping_array_1d():
 
-    F = IdentityMapping('F', dim=1)
-    f = F.get_callable_mapping()
+    f = IdentityMapping('F', dim=1)
 
     assert f.ldim == 1
     assert f.pdim == 1
@@ -241,15 +231,14 @@ def test_identity_mapping_array_1d():
     x1 = np.array([-0.7, 0.5, 33])
 
     assert np.array_equal(f(x1)[0], x1)
-    assert np.array_equal(f.jacobian(x1), np.ones((1, 1, 3)))
-    assert np.array_equal(f.jacobian_inv(x1),np.ones((1, 1, 3)))
-    assert np.array_equal(f.metric(x1), np.ones((1, 1, 3)))
-    assert np.array_equal(f.metric_det(x1), np.ones(3))
+    assert np.array_equal(f.jacobian_eval(x1), np.ones((1, 1, 3)))
+    assert np.array_equal(f.jacobian_inv_eval(x1),np.ones((1, 1, 3)))
+    assert np.array_equal(f.metric_eval(x1), np.ones((1, 1, 3)))
+    assert np.array_equal(f.metric_det_eval(x1), np.ones(3))
 
 def test_identity_mapping_array_2d():
 
-    F = IdentityMapping('F', dim=2)
-    f = F.get_callable_mapping()
+    f = IdentityMapping('F', dim=2)
 
     assert f.ldim == 2
     assert f.pdim == 2
@@ -265,16 +254,15 @@ def test_identity_mapping_array_2d():
 
     assert np.array_equal(r1, xx1)
     assert np.array_equal(r2, xx2)
-    assert np.array_equal(f.jacobian    (xx1, xx2), I)
-    assert np.array_equal(f.jacobian_inv(xx1, xx2), I)
-    assert np.array_equal(f.metric      (xx1, xx2), I)
-    assert np.array_equal(f.metric_det  (xx1, xx2), np.ones((3, 3)))
+    assert np.array_equal(f.jacobian_eval    (xx1, xx2), I)
+    assert np.array_equal(f.jacobian_inv_eval(xx1, xx2), I)
+    assert np.array_equal(f.metric_eval      (xx1, xx2), I)
+    assert np.array_equal(f.metric_det_eval  (xx1, xx2), np.ones((3, 3)))
 
 
 def test_identity_mapping_array_3d():
 
-    F = IdentityMapping('F', dim=3)
-    f = F.get_callable_mapping()
+    f = IdentityMapping('F', dim=3)
 
     assert f.ldim == 3
     assert f.pdim == 3
@@ -292,10 +280,10 @@ def test_identity_mapping_array_3d():
     assert np.array_equal(r1, xx1)
     assert np.array_equal(r2, xx2)
     assert np.array_equal(r3, xx3)
-    assert np.array_equal(f.jacobian    (xx1, xx2, xx3), I)
-    assert np.array_equal(f.jacobian_inv(xx1, xx2, xx3), I)
-    assert np.array_equal(f.metric      (xx1, xx2, xx3), I)
-    assert np.array_equal(f.metric_det  (xx1, xx2, xx3), np.ones((2, 2, 2)))
+    assert np.array_equal(f.jacobian_eval    (xx1, xx2, xx3), I)
+    assert np.array_equal(f.jacobian_inv_eval(xx1, xx2, xx3), I)
+    assert np.array_equal(f.metric_eval      (xx1, xx2, xx3), I)
+    assert np.array_equal(f.metric_det_eval  (xx1, xx2, xx3), np.ones((2, 2, 2)))
 
 
 #------------------------------------------------------------------------------
@@ -304,8 +292,7 @@ def test_affine_mapping_array_1d():
     # x = 1 - x1
     params = {'c1': 1, 'a11': -1}
 
-    F = AffineMapping('F', **params, dim=1)
-    f = F.get_callable_mapping()
+    f = AffineMapping('F', **params, dim=1)
 
     assert f.ldim == 1
     assert f.pdim == 1
@@ -313,10 +300,10 @@ def test_affine_mapping_array_1d():
     x1 = np.array([0, 0.5, 1])
 
     assert np.array_equal(f(x1)[0], [1, 0.5, 0])
-    assert np.array_equal(f.jacobian(x1), [[[-1, -1, -1]]])
-    assert np.array_equal(f.jacobian_inv(x1), [[[-1, -1, -1]]])
-    assert np.array_equal(f.metric(x1), [[[1, 1, 1]]])
-    assert np.array_equal(f.metric_det(x1), [1, 1, 1])
+    assert np.array_equal(f.jacobian_eval(x1), [[[-1, -1, -1]]])
+    assert np.array_equal(f.jacobian_inv_eval(x1), [[[-1, -1, -1]]])
+    assert np.array_equal(f.metric_eval(x1), [[[1, 1, 1]]])
+    assert np.array_equal(f.metric_det_eval(x1), [1, 1, 1])
 
 
 def test_affine_mapping_array_2d():
@@ -326,8 +313,7 @@ def test_affine_mapping_array_2d():
 
     params = dict(c1=c1, c2=c2, a11=a11, a12=a12, a21=a21, a22=a22)
 
-    F = AffineMapping('F', **params, dim=2)
-    f = F.get_callable_mapping()
+    f = AffineMapping('F', **params, dim=2)
 
     assert f.ldim == 2
     assert f.pdim == 2
@@ -357,10 +343,10 @@ def test_affine_mapping_array_2d():
 
     assert np.array_equal(r1, x)
     assert np.array_equal(r2, y)
-    assert np.array_equal(f.jacobian    (xx1, xx2), J    )
-    assert np.array_equal(f.jacobian_inv(xx1, xx2), J_inv)
-    assert np.array_equal(f.metric      (xx1, xx2), G    )
-    assert np.array_equal(f.metric_det  (xx1, xx2), G_det)
+    assert np.array_equal(f.jacobian_eval    (xx1, xx2), J    )
+    assert np.array_equal(f.jacobian_inv_eval(xx1, xx2), J_inv)
+    assert np.array_equal(f.metric_eval      (xx1, xx2), G    )
+    assert np.array_equal(f.metric_det_eval  (xx1, xx2), G_det)
 
 
 def test_affine_mapping_array_3d():
@@ -376,8 +362,7 @@ def test_affine_mapping_array_3d():
             a21=a21, a22=a22, a23=a23,
             a31=a31, a32=a32, a33=a33,)
 
-    F = AffineMapping('F', **params, dim=3)
-    f = F.get_callable_mapping()
+    f = AffineMapping('F', **params, dim=3)
 
     assert f.ldim == 3
     assert f.pdim == 3
@@ -416,10 +401,10 @@ def test_affine_mapping_array_3d():
     assert np.allclose(r1, x, rtol=RTOL, atol=ATOL)
     assert np.allclose(r2, y, rtol=RTOL, atol=ATOL)
     assert np.allclose(r3, z, rtol=RTOL, atol=ATOL)
-    assert np.array_equal(f.jacobian    (xx1, xx2, xx3), J    )
-    assert np.array_equal(f.jacobian_inv(xx1, xx2, xx3), J_inv)
-    assert np.array_equal(f.metric      (xx1, xx2, xx3), G    )
-    assert np.array_equal(f.metric_det  (xx1, xx2, xx3), G_det)
+    assert np.array_equal(f.jacobian_eval    (xx1, xx2, xx3), J    )
+    assert np.array_equal(f.jacobian_inv_eval(xx1, xx2, xx3), J_inv)
+    assert np.array_equal(f.metric_eval      (xx1, xx2, xx3), G    )
+    assert np.array_equal(f.metric_det_eval  (xx1, xx2, xx3), G_det)
 
 
 #------------------------------------------------------------------------------
@@ -430,8 +415,7 @@ def test_polar_mapping_array():
 
     params = dict(c1=c1, c2=c2, rmin=rmin, rmax=rmax)
 
-    F = PolarMapping('F', **params)
-    f = F.get_callable_mapping()
+    f = PolarMapping('F', **params)
 
     assert f.ldim == 2
     assert f.pdim == 2
@@ -463,58 +447,7 @@ def test_polar_mapping_array():
 
     assert np.array_equal(r1, x)
     assert np.array_equal(r2, y)
-    assert np.allclose(f.jacobian    (xx1, xx2), J    , rtol=RTOL, atol=ATOL)
-    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
+    assert np.allclose(f.jacobian_eval    (xx1, xx2), J    , rtol=RTOL, atol=ATOL)
+    assert np.allclose(f.jacobian_inv_eval(xx1, xx2), J_inv, rtol=RTOL, atol=ATOL)
+    assert np.allclose(f.metric_eval      (xx1, xx2), G    , rtol=RTOL, atol=ATOL)
+    assert np.allclose(f.metric_det_eval  (xx1, xx2), G_det, rtol=RTOL, atol=ATOL)
diff --git a/sympde/topology/tests/test_symbolic_expr.py b/sympde/topology/tests/test_symbolic_expr.py
index 43d1151c..cf4a7322 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 BaseMapping, 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 = BaseMapping('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 adf2cd7b..a380f5f8 100644
--- a/sympde/topology/tests/test_topology.py
+++ b/sympde/topology/tests/test_topology.py
@@ -6,7 +6,7 @@
 from sympde.topology import Boundary, NormalVector, TangentVector
 from sympde.topology import Connectivity, Edge
 from sympde.topology import Domain, ElementDomain
-from sympde.topology import Area, Mapping
+from sympde.topology import Area, BaseMapping
 from sympde.topology import Interface
 from sympde.topology import Line, Square, Cube
 from sympde.topology import IdentityMapping
@@ -43,8 +43,8 @@ def test_topology_1():
     B = Square('B')
     # ...
 
-    M1 = Mapping('M1', dim=2)
-    M2 = Mapping('M2', dim=2)
+    M1 = BaseMapping('M1', dim=2)
+    M2 = BaseMapping('M2', dim=2)
 
     D1 = M1(A)
     D2 = M2(B)
diff --git a/sympde/utilities/tests/test_utilities.py b/sympde/utilities/tests/test_utilities.py
index f6fbfe55..b2ddc9c1 100644
--- a/sympde/utilities/tests/test_utilities.py
+++ b/sympde/utilities/tests/test_utilities.py
@@ -3,7 +3,7 @@
 
 from sympy import Matrix, symbols, Array
 from sympy import S
-from sympde.utilities.utils import lambdify_sympde
+from sympde.topology.base_analytic_mapping import lambdify_sympde
 
 
 
diff --git a/sympde/utilities/utils.py b/sympde/utilities/utils.py
index 3496a1df..895ef050 100644
--- a/sympde/utilities/utils.py
+++ b/sympde/utilities/utils.py
@@ -5,86 +5,9 @@
 from mpl_toolkits.mplot3d import *
 import matplotlib.pyplot as plt
 
-from sympde.topology import IdentityMapping, InteriorDomain, MultiPatchMapping
-from sympde.topology.analytical_mapping import TorusMapping
+from sympde.topology import InteriorDomain, MultiPatchMapping
+from sympde.topology.analytic_mappings import TorusMapping, IdentityMapping
 
-def lambdify_sympde(variables, expr):
-    """
-    Custom lambify function that covers the
-    shortcomings of sympy's lambdify. Most notably,
-    this function uses numpy broadcasting rules to
-    compute the shape of the output.
-
-    Parameters
-    ----------
-    variables : sympy.core.symbol.Symbol or list of sympy.core.symbol.Symbol
-        variables that appear in the expression
-    expr :
-        Sympy expression
-
-    Returns
-    -------
-    lambda_f : callable
-        Lambdified function built using numpy.
-
-    Notes
-    -----
-    Compared to Sympy's lambdify, this function
-    is capable of properly handling constant values,
-    and array_like structures where not all components
-    depend on all variables. See below.
-
-    Examples
-    --------
-    >>> import numpy as np
-    >>> from sympy import symbols,  Matrix
-    >>> from sympde.utilities.utils import lambdify_sympde
-    >>> x, y = symbols("x,y")
-    >>> expr = Matrix([[x, x + y], [0, y]])
-    >>> f = lambdify_sympde([x,y], expr)
-    >>> f(np.array([[0, 1]]), np.array([[2], [3]]))
-    array([[[[0., 1.],
-             [0., 1.]],
-
-            [[2., 3.],
-             [3., 4.]]],
-
-
-           [[[0., 0.],
-             [0., 0.]],
-
-            [[2., 2.],
-             [3., 3.]]]])
-    """
-    array_expr = np.asarray(expr)
-    scalar_shape = array_expr.shape
-    if scalar_shape == ():
-        f = lambdify(variables, expr, 'numpy')
-        def f_vec_sc(*XYZ):
-            b = np.broadcast(*XYZ)
-            if b.ndim == 0:
-                return f(*XYZ)
-            temp = np.asarray(f(*XYZ))
-            if b.shape == temp.shape:
-                return temp
-
-            result = np.zeros(b.shape)
-            result[...] = temp
-            return result
-        return f_vec_sc
-
-    else:
-        scalar_functions = {}
-        for multi_index in it.product(*tuple(range(s) for s in scalar_shape)):
-            scalar_functions[multi_index] = lambdify(variables, array_expr[multi_index], 'numpy')
-
-        def f_vec_v(*XYZ):
-            b = np.broadcast(*XYZ)
-            result = np.zeros(scalar_shape + b.shape)
-            for multi_index in it.product(*tuple(range(s) for s in scalar_shape)):
-                result[multi_index] = scalar_functions[multi_index](*XYZ)
-            return result
-        return f_vec_v
 
 
 def plot_domain(domain, draw=True, isolines=False, refinement=None):