diff --git a/psydac/api/ast/evaluation.py b/psydac/api/ast/evaluation.py
index e4f9aa931..e70012448 100644
--- a/psydac/api/ast/evaluation.py
+++ b/psydac/api/ast/evaluation.py
@@ -1,78 +1,30 @@
-# TODO: - replace call to print_expression by SymbolicExpr (may need to use LogicalExpr)
-
 from collections import OrderedDict
-from itertools import groupby
-import numpy as np
-
-from sympy import Basic
-from sympy import symbols, Symbol, IndexedBase, Indexed, Function
-from sympy import Mul, Add, Tuple, Min, Max, Pow
-from sympy import Matrix, ImmutableDenseMatrix
-from sympy import sqrt as sympy_sqrt
-from sympy import S as sympy_S
-from sympy import Integer, Float
-from sympy.core.relational    import Le, Ge
-from sympy.logic.boolalg      import And
-from sympy import Mod, Abs
-from sympy.core.function import AppliedUndef
-
-from pyccel.ast.core import Variable, IndexedVariable
-from pyccel.ast.core import For
-from pyccel.ast.core import Assign
-from pyccel.ast.core import AugAssign
-from pyccel.ast.core import Slice
-from pyccel.ast.core import Range, Product
-from pyccel.ast.core import FunctionDef
-from pyccel.ast.core import FunctionCall
-from pyccel.ast.core import Import
-from pyccel.ast import Zeros
-from pyccel.ast import Import
-from pyccel.ast import DottedName
-from pyccel.ast import Nil
-from pyccel.ast import Len
-from pyccel.ast import If, Is, Return
-from pyccel.ast import String, Print, Shape
-from pyccel.ast import Comment, NewLine
-from pyccel.ast.core      import _atomic
-from pyccel.ast.utilities import build_types_decorator
 
+from sympy import symbols
+from sympy import Tuple, Pow
 
-from sympde.core import Cross_3d
-from sympde.core import Constant
-from sympde.calculus import grad
-from sympde.topology import Mapping
-from sympde.topology import ScalarField
-from sympde.topology import VectorField, IndexedVectorField
-from sympde.topology import Boundary, BoundaryVector, NormalVector, TangentVector
-from sympde.topology import Covariant, Contravariant
-from sympde.topology import ElementArea
-from sympde.topology import LogicalExpr
-from sympde.topology import SymbolicExpr
-from sympde.topology.derivatives import _partial_derivatives
+from sympde.topology             import Mapping
+from sympde.topology             import ScalarField
+from sympde.topology             import VectorField
+from sympde.topology             import SymbolicExpr
+from sympde.topology.space       import ScalarTestFunction
 from sympde.topology.derivatives import _logical_partial_derivatives
-from sympde.topology.derivatives import get_max_partial_derivatives
-from sympde.topology.space import ScalarFunctionSpace
-from sympde.topology.space import ScalarTestFunction
-from sympde.topology.space import VectorTestFunction
-from sympde.topology.space import IndexedTestTrial
-from sympde.topology.space import Trace
-from sympde.topology.derivatives import print_expression
-from sympde.topology.derivatives import get_atom_derivatives
-from sympde.topology.derivatives import get_index_derivatives
-from sympde.expr import BilinearForm, LinearForm, Functional, BasicForm
-
-from psydac.fem.splines import SplineSpace
-from psydac.fem.tensor  import TensorFemSpace
-from psydac.fem.vector  import ProductFemSpace
-
-from .basic import SplBasic
+
+from pyccel.ast.core      import IndexedVariable
+from pyccel.ast.core      import For
+from pyccel.ast.core      import Assign
+from pyccel.ast.core      import Slice
+from pyccel.ast.core      import Range
+from pyccel.ast.core      import FunctionDef
+from pyccel.ast.utilities import build_types_decorator
+
+from .basic     import SplBasic
 from .utilities import build_pythran_types_header, variables
 from .utilities import filter_loops, filter_product
 from .utilities import rationalize_eval_mapping
 from .utilities import compute_atoms_expr_mapping
 from .utilities import compute_atoms_expr_field
 
-
 #==============================================================================
 class EvalQuadratureMapping(SplBasic):
 
@@ -100,7 +52,7 @@ def __new__(cls, space, mapping, discrete_boundary=None, name=None,
         # ...
 
         # ...
-        ops = _partial_derivatives[:dim]
+        ops = _logical_partial_derivatives[:dim]
         M = mapping
 
         components = [M[i] for i in range(0, dim)]
@@ -196,10 +148,8 @@ def _initialize(self):
         space = self.space
         dim = space.ldim
 
-        _print = lambda i: print_expression(i, mapping_name=False)
-        mapping_atoms = [_print(f) for f in self.components]
-        mapping_str = [_print(f) for f in self.elements]
-#        mapping_str = sorted([_print(f) for f in self.elements])
+        mapping_atoms = [SymbolicExpr(f).name for f in self.components]
+        mapping_str   = [SymbolicExpr(f).name for f in self.elements  ]
 
         # ... declarations
         degrees        = variables( 'p1:%s'%(dim+1), 'int')
@@ -229,7 +179,7 @@ def _initialize(self):
 
             # ...
             nderiv = self.nderiv
-            ops = _partial_derivatives[:dim]
+            ops = _logical_partial_derivatives[:dim]
 
             if nderiv > 0:
                 weights_elements += [d(weights_pts) for d in ops]
@@ -242,7 +192,7 @@ def _initialize(self):
                 raise NotImplementedError('TODO')
             # ...
 
-            mapping_weights_str = [_print(f) for f in weights_elements]
+            mapping_weights_str = [SymbolicExpr(f).name for f in weights_elements]
             mapping_wvalues = variables(['{}_values'.format(f) for f in mapping_weights_str],
                                                 dtype='real', rank=dim, cls=IndexedVariable)
 
@@ -267,11 +217,12 @@ def _initialize(self):
         # ...
 
         # ...
-        Nj = ScalarTestFunction(space, name='Nj')
-        body = []
-        init_basis = OrderedDict()
-        updates = []
+        Nj           = ScalarTestFunction(space, name='Nj')
+        body         = []
+        init_basis   = OrderedDict()
+        updates      = []
         atomic_exprs = self.elements + weights_elements
+
         inits, updates = compute_atoms_expr_mapping(atomic_exprs, indices_quad,
                                                     indices_basis, basis, Nj)
 
@@ -356,7 +307,6 @@ def __new__(cls, space, fields, discrete_boundary=None, name=None,
         obj._backend = backend
         obj._func = obj._initialize()
 
-
         return obj
 
     @property
@@ -431,8 +381,8 @@ def _initialize(self):
 
         self._fields = Tuple(*fields)
 
-        fields_str    = [print_expression(f) for f in self._fields]
-        fields_val    = variables(['{}_values'.format(f) for f in fields_str],
+        fields_str = [SymbolicExpr(f).name for f in self._fields]
+        fields_val = variables(['{}_values'.format(f) for f in fields_str],
                                   dtype='real', rank=dim, cls=IndexedVariable)
 
 
@@ -539,7 +489,7 @@ def _initialize(self):
         indices_quad  = variables('g1:%s'%(dim+1),  'int')
         basis         = variables('basis1:%s'%(dim+1),
                                   dtype='real', rank=3, cls=IndexedVariable)
-        coeffs = ['coeff_{}'.format(print_expression(f)) for f in vector_field_atoms]
+        coeffs = ['coeff_{}'.format(SymbolicExpr(f).name) for f in vector_field_atoms]
         vector_fields_coeffs = variables(coeffs, dtype='real', rank=dim, cls=IndexedVariable)
         # ...
 
@@ -568,7 +518,7 @@ def _initialize(self):
             init_map[str(stmt.lhs)] = stmt
 
         self._vector_fields = Tuple(*vector_fields)
-        vector_fields_str   =  [print_expression(f) for f in self.vector_fields]
+        vector_fields_str   = [SymbolicExpr(f).name for f in self.vector_fields]
         vector_fields_val   = variables(['{}_values'.format(f) for f in vector_fields_str],
                                           dtype='real', rank=dim, cls=IndexedVariable)
 
@@ -626,6 +576,7 @@ def _create_loop(indices, ranges, body):
     return body
 
 #==============================================================================
+# NOTE: this is used in module 'psydac.api.ast.glt'
 class EvalArrayField(SplBasic):
 
     def __new__(cls, space, fields, discrete_boundary=None, name=None,
@@ -645,7 +596,6 @@ def __new__(cls, space, fields, discrete_boundary=None, name=None,
         obj._backend = backend
         obj._func = obj._initialize()
 
-
         return obj
 
     @property
@@ -680,7 +630,7 @@ def _initialize(self):
         mapping = self.mapping
 
         field_atoms = self.fields.atoms(ScalarField)
-        fields_str = sorted([print_expression(f) for f in self.fields])
+        fields_str = sorted([SymbolicExpr(f).name for f in self.fields])
 
         # ... declarations
         degrees        = variables( 'p1:%s'%(dim+1), 'int')
@@ -771,6 +721,7 @@ def _initialize(self):
 
 
 #==============================================================================
+# NOTE: this is used in module 'psydac.api.ast.glt'
 class EvalArrayMapping(SplBasic):
 
     def __new__(cls, space, mapping, name=None,
@@ -795,7 +746,7 @@ def __new__(cls, space, mapping, name=None,
         # ...
 
         # ...
-        ops = _partial_derivatives[:dim]
+        ops = _logical_partial_derivatives[:dim]
         M = mapping
 
         components = [M[i] for i in range(0, dim)]
@@ -880,10 +831,8 @@ def _initialize(self):
         space = self.space
         dim = space.ldim
 
-        _print = lambda i: print_expression(i, mapping_name=False)
-        mapping_atoms = [_print(f) for f in self.components]
-        mapping_str = [_print(f) for f in self.elements]
-#        mapping_str = sorted([_print(f) for f in self.elements])
+        mapping_atoms = [SymbolicExpr(f).name for f in self.components]
+        mapping_str   = [SymbolicExpr(f).name for f in self.elements  ]
 
         # ... declarations
         degrees        = variables( 'p1:%s'%(dim+1), 'int')
@@ -917,7 +866,7 @@ def _initialize(self):
 
             # ...
             nderiv = self.nderiv
-            ops = _partial_derivatives[:dim]
+            ops = _logical_partial_derivatives[:dim]
 
             if nderiv > 0:
                 weights_elements += [d(weights_pts) for d in ops]
@@ -930,7 +879,7 @@ def _initialize(self):
                 raise NotImplementedError('TODO')
             # ...
 
-            mapping_weights_str = [_print(f) for f in weights_elements]
+            mapping_weights_str = [SymbolicExpr(f).name for f in weights_elements]
             mapping_wvalues = variables(['{}_values'.format(f) for f in mapping_weights_str],
                                                 dtype='real', rank=dim, cls=IndexedVariable)
 
@@ -983,8 +932,8 @@ def _initialize(self):
         assign_spans = []
         for x, i_span, span in zip(indices_quad, i_spans, spans):
             assign_spans += [Assign(i_span, span[x])]
-
         body = assign_spans + body
+
         # put the body in for loops of quadrature points
         body = _create_loop(indices_quad, ranges_quad, body)
 
diff --git a/psydac/api/ast/expr.py b/psydac/api/ast/expr.py
index a8f2f3187..a9b9f5d91 100644
--- a/psydac/api/ast/expr.py
+++ b/psydac/api/ast/expr.py
@@ -1,20 +1,9 @@
-
-from collections import OrderedDict
-from itertools import groupby
-import string
-import random
-import numpy as np
-
-from sympy import Basic
-from sympy import symbols, Symbol, IndexedBase, Indexed, Function
-from sympy import Mul, Add, Tuple
+from sympy import symbols, Symbol, IndexedBase
+from sympy import Mul, Tuple
 from sympy import Matrix, ImmutableDenseMatrix
-from sympy import sqrt as sympy_sqrt
-from sympy import S as sympy_S
-from sympy import simplify, expand
 from sympy.core.numbers import ImaginaryUnit
 
-from pyccel.ast.core import Variable, IndexedVariable
+from pyccel.ast.core import IndexedVariable
 from pyccel.ast.core import For
 from pyccel.ast.core import Assign
 from pyccel.ast.core import AugAssign
@@ -22,42 +11,31 @@
 from pyccel.ast.core import Range
 from pyccel.ast.core import FunctionDef
 from pyccel.ast.core import FunctionCall
-from pyccel.ast import Zeros
-from pyccel.ast import Import
-from pyccel.ast import DottedName
-from pyccel.ast import Nil
-from pyccel.ast import Len
-from pyccel.ast import If, Is, Return
-from pyccel.ast import String, Print, Shape
-from pyccel.ast import Comment, NewLine
-from pyccel.ast.core      import _atomic
-
-from sympde.core import Constant
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import IndexedVectorField
-from sympde.topology import Mapping
-from sympde.expr import BilinearForm
+from pyccel.ast      import Import
+from pyccel.ast      import Nil
+from pyccel.ast      import Len
+from pyccel.ast      import If, Is, Return
+from pyccel.ast.core import _atomic
+
+from sympde.core                 import Constant
+from sympde.topology             import ScalarField, VectorField
+from sympde.topology             import IndexedVectorField
 from sympde.topology.derivatives import _partial_derivatives
 from sympde.topology.derivatives import _logical_partial_derivatives
 from sympde.topology.derivatives import get_max_partial_derivatives
-from sympde.topology.derivatives import print_expression
 from sympde.topology.derivatives import get_atom_derivatives
 from sympde.topology.derivatives import get_index_derivatives
-from sympde.topology import LogicalExpr
-from sympde.topology import SymbolicExpr
-from sympde.topology import SymbolicDeterminant
-
-from .basic import SplBasic
-from .utilities import random_string
-from .utilities import build_pythran_types_header, variables
-from .utilities import is_scalar_field, is_vector_field, is_mapping
-from .utilities import math_atoms_as_str
-#from .evaluation import EvalArrayVectorField
-from .evaluation import EvalArrayMapping, EvalArrayField
-
-from psydac.fem.splines import SplineSpace
-from psydac.fem.tensor  import TensorFemSpace
-from psydac.fem.vector  import ProductFemSpace, VectorFemSpace
+from sympde.topology             import LogicalExpr
+from sympde.topology             import SymbolicExpr
+from sympde.topology             import SymbolicDeterminant
+
+from .basic      import SplBasic
+from .utilities  import random_string
+from .utilities  import build_pythran_types_header, variables
+from .utilities  import is_scalar_field, is_vector_field
+from .utilities  import math_atoms_as_str
+
+from psydac.fem.vector import ProductFemSpace
 
 #==============================================================================
 def compute_atoms_expr(atom, basis, indices, loc_indices, dim):
@@ -268,8 +246,8 @@ def _initialize(self):
         self._fields = tuple(expr.atoms(ScalarField))
         self._vector_fields = tuple(expr.atoms(VectorField))
         # ...
-        fields_str        = tuple(map(print_expression, atomic_expr_field))
-        vector_fields_str = tuple(map(print_expression, atomic_expr_vector_field))  
+        fields_str        = tuple(SymbolicExpr(f).name for f in atomic_expr_field)
+        vector_fields_str = tuple(SymbolicExpr(f).name for f in atomic_expr_vector_field)
         
         fields = symbols(fields_str)
         vector_fields = symbols(vector_fields_str)
@@ -362,10 +340,12 @@ def _initialize(self):
                 val = val[indices]
 
                 if isinstance(expr, (Matrix, ImmutableDenseMatrix)):
-                    body += [Assign(val, expr[i_row,i_col])]
+                    rhs   = SymbolicExpr(expr[i_row,i_col])
+                    body += [Assign(val, rhs)]
 
                 else:
-                    body += [Assign(val, expr)]
+                    rhs   = SymbolicExpr(expr)
+                    body += [Assign(val, rhs)]
 
         for i in range(dim-1,-1,-1):
             x = indices[i]
diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py
index a3f7cd5ac..56540ac5a 100644
--- a/psydac/api/ast/fem.py
+++ b/psydac/api/ast/fem.py
@@ -1,18 +1,10 @@
-# TODO: - replace call to print_expression by SymbolicExpr (may need to use LogicalExpr)
-
 from collections import OrderedDict
 from itertools import groupby
 import numpy as np
 
-from sympy import Basic
-from sympy import symbols, Symbol, IndexedBase, Indexed, Function
-from sympy import Mul, Add, Tuple, Min, Max, Pow
+from sympy import symbols, Symbol, IndexedBase
+from sympy import Mul, Tuple
 from sympy import Matrix, ImmutableDenseMatrix
-from sympy import sqrt as sympy_sqrt
-from sympy import S as sympy_S
-from sympy import Integer, Float
-from sympy.core.relational    import Le, Ge
-from sympy.logic.boolalg      import And
 from sympy import Mod, Abs
 from sympy.core.function import AppliedUndef
 
@@ -24,7 +16,6 @@
 from pyccel.ast.core import Range, Product
 from pyccel.ast.core import FunctionDef
 from pyccel.ast.core import FunctionCall
-from pyccel.ast.core import Import
 from pyccel.ast import Zeros
 from pyccel.ast import Import
 from pyccel.ast import DottedName
@@ -32,49 +23,41 @@
 from pyccel.ast import Len
 from pyccel.ast import If, Is, Return
 from pyccel.ast import String, Print, Shape
-from pyccel.ast import Comment, NewLine
+from pyccel.ast import Comment
 from pyccel.ast.core      import _atomic
 from pyccel.ast.utilities import build_types_decorator
 
-from sympde.core import Cross_3d
-from sympde.core import Constant
-from sympde.calculus import grad
-from sympde.topology import Mapping
-from sympde.topology import ScalarField
-from sympde.topology import VectorField, IndexedVectorField
-from sympde.topology import Boundary, BoundaryVector, NormalVector, TangentVector
-from sympde.topology import Covariant, Contravariant
-from sympde.topology import ElementArea
-from sympde.topology import LogicalExpr
-from sympde.topology import SymbolicExpr
-from sympde.topology import UndefinedSpaceType
+from sympde.core                 import Constant
+from sympde.topology             import ScalarField
+from sympde.topology             import VectorField, IndexedVectorField
+from sympde.topology             import Boundary, BoundaryVector, NormalVector, TangentVector
+from sympde.topology             import ElementArea
+from sympde.topology             import LogicalExpr
+from sympde.topology             import SymbolicExpr
+from sympde.topology             import UndefinedSpaceType
+from sympde.topology.space       import ScalarFunctionSpace, VectorFunctionSpace
+from sympde.topology.space       import ProductSpace
+from sympde.topology.space       import ScalarTestFunction
+from sympde.topology.space       import VectorTestFunction
+from sympde.topology.space       import IndexedTestTrial
 from sympde.topology.derivatives import _partial_derivatives
 from sympde.topology.derivatives import _logical_partial_derivatives
 from sympde.topology.derivatives import get_max_partial_derivatives
-from sympde.topology.space import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology.space import ProductSpace
-from sympde.topology.space import ScalarTestFunction
-from sympde.topology.space import VectorTestFunction
-from sympde.topology.space import IndexedTestTrial
-from sympde.topology.space import Trace
-from sympde.topology.derivatives import print_expression
-from sympde.topology.derivatives import get_atom_derivatives
-from sympde.topology.derivatives import get_index_derivatives
-from sympde.expr import BilinearForm, LinearForm, Functional, BasicForm
+from sympde.expr                 import BilinearForm, LinearForm, Functional
 
 from psydac.fem.splines import SplineSpace
 from psydac.fem.tensor  import TensorFemSpace
 from psydac.fem.vector  import ProductFemSpace
 
-from .basic import SplBasic
+from .basic      import SplBasic
 from .evaluation import EvalQuadratureMapping, EvalQuadratureField, EvalQuadratureVectorField
-from .utilities import random_string
-from .utilities import build_pythran_types_header, variables
-from .utilities import compute_normal_vector, compute_tangent_vector
-from .utilities import select_loops, filter_product
-from .utilities import compute_atoms_expr
-from .utilities import is_scalar_field, is_vector_field
-from .utilities import math_atoms_as_str
+from .utilities  import random_string
+from .utilities  import build_pythran_types_header, variables
+from .utilities  import compute_normal_vector, compute_tangent_vector
+from .utilities  import select_loops, filter_product
+from .utilities  import compute_atoms_expr
+from .utilities  import is_scalar_field, is_vector_field
+from .utilities  import math_atoms_as_str
 
 
 FunctionalForms = (BilinearForm, LinearForm, Functional)
@@ -198,45 +181,26 @@ def init_loop_support(indices_elm, n_elements,
 # TODO take exponent to 1/dim
 def area_eval_mapping(mapping, area, dim, indices_quad, weight):
 
-    _print = lambda i: print_expression(i, mapping_name=False)
+    stmts = []
 
-    M = mapping
-    ops = _partial_derivatives[:dim]
+    # mapping components and their derivatives
+    ops      = _logical_partial_derivatives[:dim]
+    elements = [d(mapping[i]) for d in ops for i in range(0, dim)]
 
-    # ... mapping components and their derivatives
-    elements = [d(M[i]) for d in ops for i in range(0, dim)]
-    # ...
-
-    stmts = []
     # declarations
     stmts += [Comment('declarations')]
-    for atom in elements:
-        atom_name = _print(atom)
-        val_name = atom_name + '_values'
-        val  = IndexedBase(val_name)[indices_quad]
-
-        stmt = Assign(atom_name, val)
-        stmts += [stmt]
-
-    # ... inv jacobian
-    jac = mapping.det_jacobian
-    rdim = mapping.rdim
-    ops = _partial_derivatives[:rdim]
-    elements = [d(mapping[i]) for d in ops for i in range(0, rdim)]
     for e in elements:
-        new = print_expression(e, mapping_name=False)
-        new = Symbol(new)
-        jac = jac.subs(e, new)
-    # ...
+        lhs      = SymbolicExpr(e)
+        rhs_name = lhs.name + '_values'
+        rhs      = IndexedBase(rhs_name)[indices_quad]
+        stmts   += [Assign(lhs, rhs)]
 
-    # ...
+    # jacobian determinant
+    jac    = SymbolicExpr(mapping.det_jacobian)
     stmts += [AugAssign(area, '+', Abs(jac) * weight)]
-    # ...
 
     return stmts
 
-
-
 #==============================================================================
 # target is used when there are multiple expression (domain/boundaries)
 class Kernel(SplBasic):
@@ -567,21 +531,19 @@ def _initialize(self):
 
         # update dependencies
         self._dependencies += self.eval_fields
-                # TODO use print_expression
-
-
-        fields_str         = tuple(map(print_expression, fields))
-        fields_logical_str = [print_expression(f, logical=True) for f in
-                                     fields]
-        # ...
 
+        # TODO: remove these?
+        d_subs = dict(zip(_partial_derivatives, _logical_partial_derivatives))
+        fields_logical     = tuple(f.subs(d_subs) for f in fields)
+        fields_str         = tuple(SymbolicExpr(f).name for f in fields)
+        fields_logical_str = tuple(SymbolicExpr(f).name for f in fields_logical)
         # ...
 
         vector_field_atoms   = tuple(expr.atoms(VectorField))
         vector_fields        = []
+
         # ... create EvalQuadratureVectorField
         self._eval_vector_fields = []
-
         if atomic_expr_vector_field:
             keyfunc = lambda F: F.space.name
             data = sorted(vector_field_atoms, key=keyfunc)
@@ -606,9 +568,11 @@ def _initialize(self):
 
         # update dependencies
         self._dependencies  += self.eval_vector_fields
-        vector_fields_str    = tuple(print_expression(i) for i in  vector_fields)
-        vector_fields_logical_str = [print_expression(f, logical=True) for f in
-                                            vector_fields]
+
+        # TODO: remove these?
+        vector_fields_logical     = tuple(f.subs(d_subs) for f in vector_fields)
+        vector_fields_str         = tuple(SymbolicExpr(f).name for f in vector_fields)
+        vector_fields_logical_str = tuple(SymbolicExpr(f).name for f in vector_fields_logical)
         # ...
 
         # ... TODO add it as a method to basic class
@@ -714,7 +678,7 @@ def _initialize(self):
 
         fields_coeffs = variables(['coeff_{}'.format(f) for f in field_atoms],
                                           dtype='real', rank=dim, cls=IndexedVariable)
-        fields_val    = variables(['{}_values'.format(f) for f in fields_str],
+        fields_val    = variables(['{}_values'.format(f) for f in fields_logical_str],
                                           dtype='real', rank=dim, cls=IndexedVariable)
         fields_tmp_coeffs = variables(['tmp_coeff_{}'.format(f) for f in field_atoms],
                                               dtype='real', rank=dim, cls=IndexedVariable)
@@ -723,7 +687,7 @@ def _initialize(self):
         vector_fields_logical = symbols(vector_fields_logical_str)
 
         vector_field_atoms = [f[i] for f in vector_field_atoms for i in range(0, dim)]
-        coeffs = ['coeff_{}'.format(print_expression(f)) for f in vector_field_atoms]
+        coeffs = ['coeff_{}'.format(SymbolicExpr(f).name) for f in vector_field_atoms]
         vector_fields_coeffs = variables(coeffs, dtype='real', rank=dim, cls=IndexedVariable)
 
         vector_fields_val    = variables(['{}_values'.format(f) for f in vector_fields_str],
@@ -768,21 +732,14 @@ def _initialize(self):
         # ...
 
         # ...
-        mapping_elements = ()
-        mapping_coeffs = ()
-        mapping_values = ()
         if mapping:
-            _eval = self.eval_mapping
-            _print = lambda i: print_expression(i, mapping_name=False)
-
-            mapping_elements = [_print(i) for i in _eval.elements]
-            mapping_elements = symbols(tuple(mapping_elements))
-
-            mapping_coeffs = [_print(i) for i in _eval.mapping_coeffs]
-            mapping_coeffs = variables(mapping_coeffs, dtype='real', rank=dim, cls=IndexedVariable)
-
-            mapping_values = [_print(i) for i in _eval.mapping_values]
-            mapping_values = variables(mapping_values, dtype='real', rank=dim, cls=IndexedVariable)
+            mapping_elements = [SymbolicExpr(i) for i in self.eval_mapping.elements]
+            mapping_coeffs   = self.eval_mapping.mapping_coeffs
+            mapping_values   = self.eval_mapping.mapping_values
+        else:
+            mapping_elements = ()
+            mapping_coeffs   = ()
+            mapping_values   = ()
         # ...
 
         # ...
@@ -966,7 +923,7 @@ def _initialize(self):
 
             for i in range(start, end):
                 if not( i in zero_terms ):
-                    e = Mul(expr[i],wvol)
+                    e = SymbolicExpr(Mul(expr[i],wvol))
 
                     body.append(AugAssign(v[i],'+', e))
             # ...
diff --git a/psydac/api/ast/glt.py b/psydac/api/ast/glt.py
index a445c2463..74da6ec5f 100644
--- a/psydac/api/ast/glt.py
+++ b/psydac/api/ast/glt.py
@@ -1,69 +1,48 @@
-
 from collections import OrderedDict
 from itertools import groupby
-import string
-import random
-import numpy as np
 
-from sympy import Basic
-from sympy import symbols, Symbol, IndexedBase, Indexed, Function
-from sympy import Mul, Add, Tuple
+from sympy import symbols, Symbol, IndexedBase
+from sympy import Tuple
 from sympy import Matrix, ImmutableDenseMatrix
-from sympy import sqrt as sympy_sqrt
-from sympy import S as sympy_S
 from sympy import simplify, expand
 from sympy.core.numbers import ImaginaryUnit
 
-from pyccel.ast.core import Variable, IndexedVariable
+from pyccel.ast.core import IndexedVariable
 from pyccel.ast.core import For
 from pyccel.ast.core import Assign
-from pyccel.ast.core import AugAssign
 from pyccel.ast.core import Slice
 from pyccel.ast.core import Range
 from pyccel.ast.core import FunctionDef
 from pyccel.ast.core import FunctionCall
-from pyccel.ast import Zeros
-from pyccel.ast import Import
-from pyccel.ast import DottedName
-from pyccel.ast import Nil
-from pyccel.ast import Len
-from pyccel.ast import If, Is, Return
-from pyccel.ast import String, Print, Shape
-from pyccel.ast import Comment, NewLine
-from pyccel.ast.core      import _atomic
-
-from sympde.core import Constant
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import IndexedVectorField
-from sympde.topology import Mapping
-from sympde.expr import BilinearForm
+from pyccel.ast      import Zeros
+from pyccel.ast      import Import
+from pyccel.ast      import DottedName
+from pyccel.ast      import Nil
+from pyccel.ast      import Len
+from pyccel.ast      import If, Is, Return
+from pyccel.ast.core import _atomic
+
+from sympde.topology             import ScalarField, VectorField
+from sympde.topology             import IndexedVectorField
 from sympde.topology.derivatives import _partial_derivatives
 from sympde.topology.derivatives import _logical_partial_derivatives
 from sympde.topology.derivatives import get_max_partial_derivatives
-from sympde.topology.derivatives import print_expression
-from sympde.topology.derivatives import get_atom_derivatives
-from sympde.topology.derivatives import get_index_derivatives
-from sympde.topology import LogicalExpr
-from sympde.topology import SymbolicExpr
-from sympde.topology import SymbolicDeterminant
-
-#from gelato.core import gelatize
+from sympde.topology             import LogicalExpr
+from sympde.topology             import SymbolicExpr
+from sympde.topology             import SymbolicDeterminant
 
-from gelato.glt import BasicGlt
 from gelato.expr import gelatize
 
-from .basic import SplBasic
-from .utilities import random_string
-from .utilities import build_pythran_types_header, variables
-from .utilities import is_scalar_field, is_vector_field, is_mapping
-from .utilities import math_atoms_as_str
-#from .evaluation import EvalArrayVectorField
+from .basic      import SplBasic
+from .utilities  import random_string
+from .utilities  import build_pythran_types_header, variables
+from .utilities  import is_scalar_field, is_vector_field, is_mapping
+from .utilities  import math_atoms_as_str
 from .evaluation import EvalArrayMapping, EvalArrayField
 
-from psydac.fem.splines import SplineSpace
-from psydac.fem.tensor  import TensorFemSpace
 from psydac.fem.vector  import ProductFemSpace
 
+#==============================================================================
 class GltKernel(SplBasic):
 
     def __new__(cls, expr, spaces, name=None, mapping=None, is_rational_mapping=None, backend=None):
@@ -300,10 +279,11 @@ def _initialize(self):
         # ...
 
         # ...
-        fields_str    = sorted(tuple(map(print_expression, atomic_expr_field)))
-        fields_logical_str = sorted([print_expression(f, logical=True) for f in
-                                     atomic_expr_field])
-        field_atoms   = tuple(expr.atoms(ScalarField))
+        d_subs = dict(zip(_partial_derivatives, _logical_partial_derivatives))
+        atomic_expr_field_logical = tuple(f.subs(d_subs) for f in atomic_expr_field)
+        fields_str         = tuple(sorted(SymbolicExpr(f).name for f in atomic_expr_field))
+        fields_logical_str = tuple(sorted(SymbolicExpr(f).name for f in atomic_expr_field_logical))
+        field_atoms        = tuple(expr.atoms(ScalarField))
         # ...
 
         # ... create EvalArrayField
@@ -381,7 +361,7 @@ def _initialize(self):
 #        vector_fields_logical = symbols(vector_fields_logical_str)
 #
 #        vector_field_atoms = [f[i] for f in vector_field_atoms for i in range(0, dim)]
-#        coeffs = ['coeff_{}'.format(print_expression(f)) for f in vector_field_atoms]
+#        coeffs = ['coeff_{}'.format(SymbolicExpr(f).name) for f in vector_field_atoms]
 #        vector_fields_coeffs = variables(coeffs, dtype='real', rank=dim, cls=IndexedVariable)
 #
 #        vector_fields_val    = variables(['{}_values'.format(f) for f in vector_fields_str],
@@ -389,21 +369,14 @@ def _initialize(self):
         # ...
 
         # ...
-        mapping_elements = ()
-        mapping_coeffs = ()
-        mapping_values = ()
         if mapping:
-            _eval = self.eval_mapping
-            _print = lambda i: print_expression(i, mapping_name=False)
-
-            mapping_elements = [_print(i) for i in _eval.elements]
-            mapping_elements = symbols(tuple(mapping_elements))
-
-            mapping_coeffs = [_print(i) for i in _eval.mapping_coeffs]
-            mapping_coeffs = variables(mapping_coeffs, dtype='real', rank=dim, cls=IndexedVariable)
-
-            mapping_values = [_print(i) for i in _eval.mapping_values]
-            mapping_values = variables(mapping_values, dtype='real', rank=dim, cls=IndexedVariable)
+            mapping_elements = [SymbolicExpr(i) for i in self.eval_mapping.elements]
+            mapping_coeffs   = self.eval_mapping.mapping_coeffs
+            mapping_values   = self.eval_mapping.mapping_values
+        else:
+            mapping_elements = ()
+            mapping_coeffs   = ()
+            mapping_values   = ()
         # ...
 
 #        # ...
@@ -457,10 +430,12 @@ def _initialize(self):
                 symbol = symbol[indices]
 
                 if isinstance(expr, (Matrix, ImmutableDenseMatrix)):
-                    body += [Assign(symbol, expr[i_row,i_col])]
+                    rhs   = SymbolicExpr(expr[i_row,i_col])
+                    body += [Assign(symbol, rhs)]
 
                 else:
-                    body += [Assign(symbol, expr)]
+                    rhs   = SymbolicExpr(expr)
+                    body += [Assign(symbol, rhs)]
 
         for i in range(dim-1,-1,-1):
             x = indices[i]
diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py
index 4f0c2f3c5..b0575685d 100644
--- a/psydac/api/ast/utilities.py
+++ b/psydac/api/ast/utilities.py
@@ -4,38 +4,36 @@
 import numpy as np
 
 from sympy import Symbol, IndexedBase, Indexed
-from sympy import Mul, Add, Tuple, Min, Max, Pow
+from sympy import Mul, Tuple, Pow
 from sympy import Matrix
 from sympy import sqrt as sympy_sqrt
 from sympy.utilities.iterables import cartes
 
+from sympde.topology.space       import ScalarTestFunction
+from sympde.topology.space       import VectorTestFunction
+from sympde.topology.space       import IndexedTestTrial
+from sympde.topology             import ScalarField
+from sympde.topology             import VectorField, IndexedVectorField
+from sympde.topology             import Mapping
+from sympde.topology.derivatives import _partial_derivatives
+from sympde.topology.derivatives import _logical_partial_derivatives
+from sympde.topology.derivatives import get_atom_derivatives
+from sympde.topology.derivatives import get_index_derivatives
+from sympde.topology.derivatives import get_atom_logical_derivatives
+from sympde.topology.derivatives import get_index_logical_derivatives
+from sympde.topology             import LogicalExpr
+from sympde.topology             import SymbolicExpr
+from sympde.core                 import Cross_3d
+
 from pyccel.ast.core import Variable, IndexedVariable
 from pyccel.ast.core import For
 from pyccel.ast.core import Assign
 from pyccel.ast.core import AugAssign
 from pyccel.ast.core import Range, Product
 from pyccel.ast.core import _atomic
-from pyccel.ast import Comment, NewLine
-
-from psydac.api.printing.pycode      import pycode
-
-
-from sympde.topology.derivatives import _partial_derivatives
-from sympde.topology.derivatives import _logical_partial_derivatives
-from sympde.topology.space import ScalarTestFunction
-from sympde.topology.space import VectorTestFunction
-from sympde.topology.space import IndexedTestTrial
-from sympde.topology import ScalarField
-from sympde.topology import VectorField, IndexedVectorField
-from sympde.topology import Mapping
-from sympde.topology.derivatives import print_expression
-from sympde.topology.derivatives import get_atom_derivatives
-from sympde.topology.derivatives import get_index_derivatives
-from sympde.topology import LogicalExpr
-from sympde.topology import SymbolicExpr
-from sympde.core import Cross_3d
-
+from pyccel.ast      import Comment
 
+from psydac.api.printing.pycode import pycode  # TODO: remove
 
 #==============================================================================
 def random_string( n ):
@@ -64,6 +62,9 @@ def is_scalar_field(expr):
     if isinstance(expr, _partial_derivatives):
         return is_scalar_field(expr.args[0])
 
+    elif isinstance(expr, _logical_partial_derivatives):
+        return is_scalar_field(expr.args[0])
+
     elif isinstance(expr, ScalarField):
         return True
 
@@ -75,6 +76,9 @@ def is_vector_field(expr):
     if isinstance(expr, _partial_derivatives):
         return is_vector_field(expr.args[0])
 
+    elif isinstance(expr, _logical_partial_derivatives):
+        return is_vector_field(expr.args[0])
+
     elif isinstance(expr, (VectorField, IndexedVectorField)):
         return True
 
@@ -150,8 +154,8 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test,
 
     Returns
     -------
-    assignments : <list>
-       list of assignments of the atomic expression evaluated in the quadature points
+    inits : <list>
+       list of assignments of the atomic expression evaluated in the quadrature points
 
     map_stmts : <list>
         list of assigments of atomic expression in case of mapping
@@ -171,35 +175,53 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test,
     for atom in atomic_exprs:   
         if not isinstance(atom, cls):
             raise TypeError('atom must be of type {}'.format(str(cls)))
-    # ...
     
-    atomic_exprs = list(atomic_exprs)
-    new_atoms    = []
-    # ... map basis function
-    map_stmts = []
+    # If there is a mapping, compute [dx(u), dy(u), dz(u)] as functions
+    # of [dx1(u), dx2(u), dx3(u)], and store results into intermediate
+    # variables [u_x, u_y, u_z]. (Same thing is done for higher derivatives.)
+    #
+    # Accordingly, we create a new list of atoms where all partial derivatives
+    # are taken with respect to the logical coordinates.
     if mapping:
+
+        new_atoms = set()
+        map_stmts = []
+        get_index = get_index_logical_derivatives
+        get_atom  = get_atom_logical_derivatives
+
         for atom in atomic_exprs:
-             if isinstance(atom, _partial_derivatives):
-                name = print_expression(atom)
-                rhs = LogicalExpr(mapping, atom)
+
+            if isinstance(atom, _partial_derivatives):
+                lhs   = SymbolicExpr(atom)
+                rhs_p = LogicalExpr(mapping, atom)
+
                 # we look for new_atoms that must be added to atomic_exprs
                 # because we need them in the maps stmts
-                logical_atoms = _atomic(rhs, cls=_logical_partial_derivatives)
-                rhs           = SymbolicExpr(rhs)
-                map_stmts     += [Assign(Symbol(name), rhs)]
-                for atom in logical_atoms:
-                    ls = _atomic(atom, Symbol)
+                logical_atoms = _atomic(rhs_p, cls=_logical_partial_derivatives)
+                for a in logical_atoms:
+                    ls = _atomic(a, Symbol)
                     assert len(ls) == 1
                     if isinstance(ls[0], cls):
-                        new_atoms += [logical2physical(atom)]
+                        new_atoms.add(a)
     
-    atomic_exprs = {*atomic_exprs, *new_atoms}
+                rhs = SymbolicExpr(rhs_p)
+                map_stmts += [Assign(lhs, rhs)]
 
-    assigns = []
-    for atom in atomic_exprs:
+            else:
+                new_atoms.add(atom)
 
-        orders = [*get_index_derivatives(atom).values()]
-        a      = get_atom_derivatives(atom)
+    else:
+        new_atoms = atomic_exprs
+        map_stmts = []
+        get_index = get_index_derivatives
+        get_atom  = get_atom_derivatives
+
+    # Create a list of statements for initialization of the point values,
+    # for each of the atoms in our (possibly new) list.
+    inits = []
+    for atom in new_atoms:
+        orders = [*get_index(atom).values()]
+        a      = get_atom(atom)
         test   = _get_name(a) in [_get_name(f) for f in test_function]
 
         if test or is_linear:
@@ -209,16 +231,14 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test,
             basis  = basis_trial
             idxs   = indices_trial
 
-        args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
-
-        # ... assign basis on quad point
-        logical  = not( mapping is None )
-        name     = print_expression(atom, logical=logical)
-        assigns += [Assign(Symbol(name), Mul(*args))]
-
-    # ...
-    return assigns, map_stmts
+        args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
+        lhs    = SymbolicExpr(atom)
+        rhs    = Mul(*args)
+        inits += [Assign(lhs, rhs)]
 
+    # Return the initialization statements, and the additional initialization
+    # of intermediate variables in case of mapping
+    return inits, map_stmts
 
 #==============================================================================
 def compute_atoms_expr_field(atomic_exprs, indices_quad,
@@ -259,10 +279,10 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad,
         list of augmented assignments which are updated in each loop iteration
 
     map_stmts : <list>
-        list of assigments of atomic expression in case of mapping
+        list of assignments of atomic expression in case of mapping
 
-    atomic_exprs: <list>
-        list of new atomic expressions that were introduced in case of a mapping
+    new_atoms: <list>
+        updated list of atomic expressions (some were introduced in case of a mapping)
     """
 
     inits     = []
@@ -274,68 +294,84 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad,
            IndexedVectorField,
            VectorField)
 
+    # If there is a mapping, compute [dx(u), dy(u), dz(u)] as functions
+    # of [dx1(u), dx2(u), dx3(u)], and store results into intermediate
+    # variables [u_x, u_y, u_z]. (Same thing is done for higher derivatives.)
+    #
+    # Accordingly, we create a new list of atoms where all partial derivatives
+    # are taken with respect to the logical coordinates.
     if mapping:
-        # ... map basis function
-        new_atoms = []
+
+        new_atoms = set()
+        map_stmts = []
+        get_index = get_index_logical_derivatives
+        get_atom  = get_atom_logical_derivatives
+
         for atom in atomic_exprs:
+
             if isinstance(atom, _partial_derivatives):
-                name       = print_expression(atom)
-                rhs        = LogicalExpr(mapping, atom)
+                lhs   = SymbolicExpr(atom)
+                rhs_p = LogicalExpr(mapping, atom)
 
-                logical_atoms = _atomic(rhs, cls=_logical_partial_derivatives)
-                rhs           = SymbolicExpr(rhs)
-                sym           = Symbol(name)
-                map_stmts    += [Assign(sym, rhs)]
+                # we look for new_atoms that must be added to atomic_exprs
+                # because we need them in the maps stmts
+                logical_atoms = _atomic(rhs_p, cls=_logical_partial_derivatives)
+                for a in logical_atoms:
+                    ls = _atomic(a, Symbol)
+                    assert len(ls) == 1
+                    if isinstance(ls[0], cls):
+                        new_atoms.add(a)
 
-                if not pycode(atom) == name:
-                    var        = Symbol(pycode(atom))
-                    map_stmts += [Assign(var, sym)]
+                rhs = SymbolicExpr(rhs_p)
+                map_stmts += [Assign(lhs, rhs)]
 
-                for atom in logical_atoms:
-                    ls = _atomic(atom, Symbol)
-                    if isinstance(ls[0], cls):
-                        new_atoms += [logical2physical(atom)]
+            else:
+                new_atoms.add(atom)
+
+    else:
+        new_atoms = atomic_exprs
+        map_stmts = []
+        get_index = get_index_derivatives
+        get_atom  = get_atom_derivatives
 
-        atomic_exprs = {*atomic_exprs, *new_atoms}
+    # Make sure that we only pick one between 'dx1(dx2(u))' and 'dx2(dx1(u))'
+    new_atoms = {SymbolicExpr(a).name : a for a in new_atoms}
+    new_atoms = tuple(new_atoms.values())
 
-    atomic_exprs =  {print_expression(a):a for a in atomic_exprs}
-    atomic_exprs = tuple(atomic_exprs.values())
+    # Create a list of statements for initialization of the point values,
+    # for each of the atoms in our (possibly new) list.
+    inits = []
+    for atom in new_atoms:
 
-    for atom in atomic_exprs:
+        # Extract field, compute name of coefficient variable, and get base
         if is_scalar_field(atom):
-            field      = list(atom.atoms(ScalarField))[0]
-            field_name = 'coeff_' + print_expression(field)
+            field      = atom.atoms(ScalarField).pop()
+            field_name = 'coeff_' + SymbolicExpr(field).name
             base       = field
-
         elif is_vector_field(atom):
-            field      = list(atom.atoms(IndexedVectorField))[0]
-            field_name = 'coeff_' + print_expression(field)
+            field      = atom.atoms(IndexedVectorField).pop()
+            field_name = 'coeff_' + SymbolicExpr(field).name
             base       = field.base
-
         else:
             raise TypeError('atom must be either scalar or vector field')
 
-        # ...
-
-        test_fun      = atom.subs(base, test_function)
-        name          = print_expression(test_fun)
-        test_fun      = Symbol(name)
-        # ...
+        # Obtain variable for storing point values of test function
+        test_fun = SymbolicExpr(atom.subs(base, test_function))
 
-        # ...
-        orders = [*get_index_derivatives(atom).values()]
+        # ...
+        orders = [*get_index(atom).values()]
         args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
         inits += [Assign(test_fun, Mul(*args))]
         # ...
 
         # ...
         args     = [IndexedBase(field_name)[idxs], test_fun]
-        val_name = print_expression(atom) + '_values'
+        val_name = SymbolicExpr(atom).name + '_values'
         val      = IndexedBase(val_name)[indices_quad]
         updates += [AugAssign(val,'+',Mul(*args))]
         # ...
 
-    return inits, updates, map_stmts, atomic_exprs
+    return inits, updates, map_stmts, new_atoms
 
 #==============================================================================
 # TODO: merge into 'compute_atoms_expr_field'
@@ -374,30 +410,29 @@ def compute_atoms_expr_mapping(atomic_exprs, indices_quad,
 
     """
 
-    _print  = lambda i: print_expression(i, mapping_name=False)
     inits   = []
     updates = []
     for atom in atomic_exprs:
-        element = get_atom_derivatives(atom)
-        element_name = 'coeff_' + _print(element)
+
+        element = get_atom_logical_derivatives(atom)
+        element_name = 'coeff_' + SymbolicExpr(element).name
 
         # ...
         test_fun = atom.subs(element, test_function)
-        name     = print_expression(test_fun, logical=True)
-        test_fun = Symbol(name)
+        test_fun = SymbolicExpr(test_fun)
         # ...
 
         # ...
-        orders = [*get_index_derivatives(atom).values()]
+        orders = [*get_index_logical_derivatives(atom).values()]
         args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
         inits += [Assign(test_fun, Mul(*args))]
         # ...
 
         # ...
-        args     = [IndexedBase(element_name)[idxs], test_fun]
-        val_name = _print(atom) + '_values'
+        val_name = SymbolicExpr(atom).name + '_values'
         val      = IndexedBase(val_name)[indices_quad]
-        updates += [AugAssign(val,'+',Mul(*args))]
+        expr     = IndexedBase(element_name)[idxs] * test_fun
+        updates += [AugAssign(val, '+', expr)]
         # ...
 
     return inits, updates
@@ -405,11 +440,9 @@ def compute_atoms_expr_mapping(atomic_exprs, indices_quad,
 #==============================================================================
 def rationalize_eval_mapping(mapping, nderiv, space, indices_quad):
 
-    _print = lambda i: print_expression(i, mapping_name=False)
-
     M = mapping
     dim = space.ldim
-    ops = _partial_derivatives[:dim]
+    ops = _logical_partial_derivatives[:dim]
 
     # ... mapping components and their derivatives
     components = [M[i] for i in range(0, dim)]
@@ -447,7 +480,7 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad):
     # declarations
     stmts += [Comment('declarations')]
     for atom in elements + weights_elements:
-        atom_name = _print(atom)
+        atom_name = SymbolicExpr(atom).name
         val_name = atom_name + '_values'
         val  = IndexedBase(val_name)[indices_quad]
 
@@ -458,9 +491,9 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad):
     stmts += [Comment('rationalize')]
 
     # 0 order terms
-    w  = Symbol( _print( weights ) )
     for i in range(dim):
-        u  = Symbol( _print( M[i] ) )
+        w = SymbolicExpr(weights)
+        u = SymbolicExpr(M[i])
 
         val_name = u.name + '_values'
         val  = IndexedBase(val_name)[indices_quad]
@@ -471,12 +504,12 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad):
     # 1 order terms
     if nderiv >= 1:
         for d in ops:
-            w  = Symbol( _print( weights    ) )
-            dw = Symbol( _print( d(weights) ) )
+            w  = SymbolicExpr(  weights )
+            dw = SymbolicExpr(d(weights))
 
             for i in range(dim):
-                u  = Symbol( _print( M[i]    ) )
-                du = Symbol( _print( d(M[i]) ) )
+                u  = SymbolicExpr(  M[i] )
+                du = SymbolicExpr(d(M[i]))
 
                 val_name = du.name + '_values'
                 val  = IndexedBase(val_name)[indices_quad]
@@ -488,16 +521,16 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad):
     if nderiv >= 2:
         for e, d1 in enumerate(ops):
             for d2 in ops[:e+1]:
-                w     = Symbol( _print(       weights   ) )
-                d1w   = Symbol( _print(    d1(weights)  ) )
-                d2w   = Symbol( _print(    d2(weights)  ) )
-                d1d2w = Symbol( _print( d1(d2(weights)) ) )
+                w     = SymbolicExpr(      weights  )
+                d1w   = SymbolicExpr(   d1(weights) )
+                d2w   = SymbolicExpr(   d2(weights) )
+                d1d2w = SymbolicExpr(d1(d2(weights)))
 
                 for i in range(dim):
-                    u     = Symbol( _print(       M[i]   ) )
-                    d1u   = Symbol( _print(    d1(M[i])  ) )
-                    d2u   = Symbol( _print(    d2(M[i])  ) )
-                    d1d2u = Symbol( _print( d1(d2(M[i])) ) )
+                    u     = SymbolicExpr(      M[i]  )
+                    d1u   = SymbolicExpr(   d1(M[i]) )
+                    d2u   = SymbolicExpr(   d2(M[i]) )
+                    d1d2u = SymbolicExpr(d1(d2(M[i])))
 
                     val_name = d1d2u.name + '_values'
                     val  = IndexedBase(val_name)[indices_quad]
diff --git a/psydac/api/expr.py b/psydac/api/expr.py
index a70192fc0..932e794dd 100644
--- a/psydac/api/expr.py
+++ b/psydac/api/expr.py
@@ -1,25 +1,23 @@
 # coding: utf-8
 
 # TODO for the moment we assume Product of same space
+# TODO properly treat expression with mapping
 
-import numpy as np
 from itertools import product
+from sympy import Expr
+import numpy as np
+
+from sympde.expr import TerminalExpr
 
 from psydac.api.basic         import BasicCodeGen
 from psydac.api.settings      import PSYDAC_BACKEND_PYTHON, PSYDAC_DEFAULT_FOLDER
 from psydac.api.grid          import CollocationBasisValues
-
 from psydac.api.ast.expr      import ExprKernel, ExprInterface
-
 from psydac.cad.geometry      import Geometry
 from psydac.mapping.discrete  import SplineMapping, NurbsMapping
-
-from psydac.fem.splines import SplineSpace
-from psydac.fem.tensor  import TensorFemSpace
-from psydac.fem.vector  import ProductFemSpace
-
-from sympde.expr        import TerminalExpr
-from sympy              import Expr
+from psydac.fem.splines       import SplineSpace
+from psydac.fem.tensor        import TensorFemSpace
+from psydac.fem.vector        import ProductFemSpace
 
 #==============================================================================
 class DiscreteExpr(BasicCodeGen):
diff --git a/psydac/api/printing/pycode.py b/psydac/api/printing/pycode.py
index fc102b80c..a00314d0f 100644
--- a/psydac/api/printing/pycode.py
+++ b/psydac/api/printing/pycode.py
@@ -1,14 +1,11 @@
 from sympy.core import Symbol
-from sympy import Tuple
 
 from pyccel.codegen.printing.pycode import PythonCodePrinter as PyccelPythonCodePrinter
 
-from sympde.calculus import Dot, Inner, Cross
-from sympde.calculus import Grad, Rot, Curl, Div
-from sympde.topology import Line, Square, Cube
 from sympde.topology.derivatives import _partial_derivatives
-from sympde.topology.derivatives import print_expression
+from sympde.topology             import SymbolicExpr
 
+#==============================================================================
 class PythonCodePrinter(PyccelPythonCodePrinter):
 
     def __init__(self, settings=None):
@@ -59,51 +56,6 @@ def _print_GltInterface(self, expr):
         code = '\n'.join(self._print(i) for i in expr.imports)
 
         return code +'\n' + self._print(expr.func)
-    # .........................................................
-
-    # .........................................................
-    #         SYMPDE objects
-    # .........................................................
-    def _print_dx(self, expr):
-        arg = expr.args[0]
-        if isinstance(arg, _partial_derivatives):
-            arg = print_expression(arg, mapping_name=False)
-
-        else:
-            arg = self._print(arg) + '_'
-
-        return arg + 'x'
-
-    def _print_dy(self, expr):
-        arg = expr.args[0]
-        if isinstance(arg, _partial_derivatives):
-            arg = print_expression(arg, mapping_name=False)
-
-        else:
-            arg = self._print(arg) + '_'
-
-        return arg + 'y'
-
-    def _print_dz(self, expr):
-        arg = expr.args[0]
-        if isinstance(arg, _partial_derivatives):
-            arg = print_expression(arg, mapping_name=False)
-
-        else:
-            arg = self._print(arg) + '_'
-
-        return arg + 'z'
-
-    def _print_IndexedTestTrial(self, expr):
-        base = self._print(expr.base)
-        index = self._print(expr.indices[0])
-        return  '{base}_{i}'.format(base=base, i=index)
-
-    def _print_IndexedVectorField(self, expr):
-        base = self._print(expr.base)
-        index = self._print(expr.indices[0])
-        return  '{base}_{i}'.format(base=base, i=index)
-    # .........................................................
 
     # .........................................................
     #        SYMPY objects
@@ -115,9 +67,8 @@ def _print_AppliedUndef(self, expr):
         args = ','.join(self._print(i) for i in expr.args)
         fname = self._print(expr.func.__name__)
         return '{fname}({args})'.format(fname=fname, args=args)
-    # .........................................................
-
 
+#==============================================================================
 def pycode(expr, **settings):
     """ Converts an expr to a string of Python code
     Parameters
diff --git a/psydac/api/tests/test_api_1d_compatible_spaces.py b/psydac/api/tests/test_api_1d_compatible_spaces.py
index 14d1f62fe..9c8d83ab2 100644
--- a/psydac/api/tests/test_api_1d_compatible_spaces.py
+++ b/psydac/api/tests/test_api_1d_compatible_spaces.py
@@ -24,7 +24,6 @@
 from psydac.api.discretization import discretize
 
 from numpy import linspace, zeros, allclose
-import numpy as np
 from mpi4py import MPI
 import pytest
 
@@ -97,9 +96,7 @@ def run_system_1_1d_dir(f0, sol, ncells, degree):
 def test_api_system_1_1d_dir_1():
 
     from sympy.abc import x
-    
+
     f0 = -(2*pi)**2*sin(2*pi*x)
     u  = sin(2*pi*x)
     x  = run_system_1_1d_dir(f0, u,ncells=[10], degree=[2])
-
-    
diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py
index 3bd104aaf..6c73888b5 100644
--- a/psydac/api/tests/test_api_2d_scalar_mapping.py
+++ b/psydac/api/tests/test_api_2d_scalar_mapping.py
@@ -264,6 +264,72 @@ def run_laplace_2d_neu(filename, solution, f, comm=None):
 
     return l2_error, h1_error
 
+#==============================================================================
+def run_biharmonic_2d_dir(filename, solution, f, comm=None):
+
+    # ... abstract model
+    domain = Domain.from_file(filename)
+
+    V = ScalarFunctionSpace('V', domain)
+
+    F = element_of(V, name='F')
+
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
+
+    int_0 = lambda expr: integral(domain , expr)
+
+    expr = laplace(v) * laplace(u)
+    a = BilinearForm((v,u),int_0(expr))
+
+    expr = f*v
+    l = LinearForm(v, int_0(expr))
+
+    error = F - solution
+    l2norm = Norm(error, domain, kind='l2')
+    h1norm = Norm(error, domain, kind='h1')
+    h2norm = Norm(error, domain, kind='h2')
+
+    nn = NormalVector('nn')
+    bc  = [EssentialBC(u, 0, domain.boundary)]
+    bc += [EssentialBC(dot(grad(u), nn), 0, domain.boundary)]
+    equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc)
+    # ...
+
+    # ... create the computational domain from a topological domain
+    domain_h = discretize(domain, filename=filename, comm=comm)
+    # ...
+
+    # ... discrete spaces
+    Vh = discretize(V, domain_h)
+    # ...
+
+    # ... dsicretize the equation using Dirichlet bc
+    equation_h = discretize(equation, domain_h, [Vh, Vh])
+    # ...
+
+    # ... discretize norms
+    l2norm_h = discretize(l2norm, domain_h, Vh)
+    h1norm_h = discretize(h1norm, domain_h, Vh)
+    h2norm_h = discretize(h2norm, domain_h, Vh)
+    # ...
+
+    # ... solve the discrete equation
+    x = equation_h.solve()
+    # ...
+
+    # ...
+    phi = FemField( Vh, x )
+    # ...
+
+    # ... compute norms
+    l2_error = l2norm_h.assemble(F=phi)
+    h1_error = h1norm_h.assemble(F=phi)
+    h2_error = h2norm_h.assemble(F=phi)
+    # ...
+
+    return l2_error, h1_error, h2_error
+
 ###############################################################################
 #            SERIAL TESTS
 ###############################################################################
@@ -521,6 +587,27 @@ def test_api_laplace_2d_neu_identity():
     assert( abs(l2_error - expected_l2_error) < 1.e-7)
     assert( abs(h1_error - expected_h1_error) < 1.e-7)
 
+#==============================================================================
+def test_api_biharmonic_2d_dir_identity():
+    filename = os.path.join(mesh_dir, 'quart_circle.h5')
+
+    from sympy.abc import x,y
+
+    solution = (sin(pi*x)*sin(pi*y))**2
+    f        = laplace(laplace(solution))
+
+    l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f)
+
+    expected_l2_error = 0.3092792236008727
+    expected_h1_error = 1.3320441589030887
+    expected_h2_error = 6.826223014197719
+
+    assert( abs(l2_error - expected_l2_error) < 1.e-7)
+    assert( abs(h1_error - expected_h1_error) < 1.e-7)
+    assert( abs(h2_error - expected_h2_error) < 1.e-7)
+
+# TODO: add biharmonic on Collela and quart_circle mappings
+
 ##==============================================================================
 ## TODO DEBUG, not working since merge with devel
 #def test_api_laplace_2d_neu_collela():
diff --git a/psydac/api/tests/test_api_expr_2d_scalar.py b/psydac/api/tests/test_api_expr_2d_scalar.py
index d6aef4024..b88da5720 100644
--- a/psydac/api/tests/test_api_expr_2d_scalar.py
+++ b/psydac/api/tests/test_api_expr_2d_scalar.py
@@ -69,6 +69,6 @@ def run_poisson_2d_dir(ncells, degree, comm=None):
 ###############################################################################
 
 #==============================================================================
-def test_api_glt_poisson_2d_dir_1():
+def test_api_expr_2d_1():
     run_poisson_2d_dir(ncells=[2**3,2**3], degree=[2,2])