diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py
index 5c8b658..55aeee8 100644
--- a/sympde/topology/mapping.py
+++ b/sympde/topology/mapping.py
@@ -14,6 +14,7 @@
 from sympy                 import sqrt, symbols
 from sympy.core.exprtools  import factor_terms
 from sympy.polys.polytools import parallel_poly_from_expr
+from sympy.core.compatibility import is_sequence
 
 from sympde.core              import Constant
 from sympde.core.basic        import BasicMapping
@@ -895,10 +896,10 @@ def eval(cls, expr, domain, **options):
         """."""
 
         from sympde.expr.evaluation import TerminalExpr, DomainExpression
-        from sympde.expr.expr import BilinearForm, LinearForm, BasicForm, Norm
-        from sympde.expr.expr import Integral
+        from sympde.expr.expr import BilinearForm, LinearForm, Norm, Integral
 
-        types = (ScalarFunction, VectorFunction, DifferentialOperator, Trace, Integral)
+        types = (ScalarFunction, VectorFunction, DifferentialOperator,
+                 Trace, Integral, DiffOperator)
 
         mapping   = domain.mapping
         dim       = domain.dim
@@ -909,20 +910,17 @@ def eval(cls, expr, domain, **options):
         ph_coords = ['x', 'y', 'z']
 
         if not has(expr, types):
-            if has(expr, DiffOperator):
-                return cls( expr, domain, evaluate=False)
-            else:
-                syms = symbols(ph_coords[:dim])
-                if isinstance(mapping, InterfaceMapping):
-                    mapping = mapping.minus
-                    # here we assume that the two mapped domains
-                    # are identical in the interface so we choose one of them
-                Ms   = [mapping[i] for i in range(dim)]
-                expr = expr.subs(list(zip(syms, Ms)))
+            syms = symbols(ph_coords[:dim])
+            if isinstance(mapping, InterfaceMapping):
+                mapping = mapping.minus
+                # here we assume that the two mapped domains
+                # are identical in the interface so we choose one of them
+            Ms   = [mapping[i] for i in range(dim)]
+            expr = expr.subs(list(zip(syms, Ms)))
 
-                if mapping.is_analytical:
-                    expr = expr.subs(list(zip(Ms, mapping.expressions)))
-                return expr
+            if mapping.is_analytical:
+                expr = expr.subs(list(zip(Ms, mapping.expressions)))
+            return expr
 
         if isinstance(expr, Symbol) and expr.name in l_coords:
             return expr