diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py
index 560024a6..0f0ff9cd 100644
--- a/sympde/topology/derivatives.py
+++ b/sympde/topology/derivatives.py
@@ -520,22 +520,31 @@ def eval(cls, *_args):
 #==============================================================================
 class GradBasic(CalculusFunction):
 
-    nargs = None
+    nargs = None  # TODO [YG, 21.07.2022]: Check why we need this
     name = 'Grad'
 
     @cacheit
-    def __new__(cls, *args, **options):
+    def __new__(cls, u, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
-            r = cls.eval(*args)
+            return cls.eval(u)
         else:
-            r = None
+            return Basic.__new__(cls, u, **options)
 
-        if r is None:
-            return Basic.__new__(cls, *args, **options)
+    @classmethod
+    @cacheit
+    def eval(cls, u):
+        du = [dxi(u) for dxi in cls._derivatives]
+
+        # 1D case: scalar output
+        if len(du) == 1: return du[0]
+
+        # 2D and 3D cases
+        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
+            lines = [list(d[:]) for d in du]
         else:
-            return r
+            lines = [[d] for d in du]
+        return ImmutableDenseMatrix(lines)
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -544,78 +553,30 @@ def __getitem__(self, indices, **kw_args):
         else:
             return Indexed(self, indices, **kw_args)
 
-class Grad_1d(GradBasic):
-
-    @classmethod
-    @cacheit
-    def eval(cls, *_args):
-
-        if not _args:
-            return
 
-        u = _args[0]
+class Grad_1d(GradBasic):
+    _derivatives = (dx,)
 
-        return dx(u)
 
 class Grad_2d(GradBasic):
+    _derivatives = (dx, dy)
 
-    @classmethod
-    @cacheit
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u  = _args[0]
-        du = (dx(u), dy(u))
-
-        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
-            lines = [list(d[:]) for d in du]
-        else:
-            lines = [[d] for d in du]
-
-        v = ImmutableDenseMatrix(lines)
-        return v
 
 class Grad_3d(GradBasic):
-
-    @classmethod
-    @cacheit
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u  = _args[0]
-        du = (dx(u), dy(u), dz(u))
-
-        if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)):
-            lines = [list(d[:]) for d in du]
-        else:
-            lines = [[d] for d in du]
-
-        v = ImmutableDenseMatrix(lines)
-
-        return v
+    _derivatives = (dx, dy, dz)
 
 #==============================================================================
 class CurlBasic(CalculusFunction):
 
-    nargs = None
+    nargs = None  # TODO [YG, 21.07.2022]: Check why we need this
     name = 'Curl'
 
-    def __new__(cls, *args, **options):
+    def __new__(cls, u, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
-            r = cls.eval(*args)
-        else:
-            r = None
-
-        if r is None:
-            return Basic.__new__(cls, *args, **options)
+            return cls.eval(u)
         else:
-            return r
+            return Basic.__new__(cls, u, **options)
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -624,28 +585,18 @@ def __getitem__(self, indices, **kw_args):
         else:
             return Indexed(self, indices, **kw_args)
 
+
 class Curl_2d(CurlBasic):
 
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
+    def eval(cls, u):
+        return dx(u[1]) - dy(u[0])
 
-        return dx(u[1])-dy(u[0])
 
 class Curl_3d(CurlBasic):
 
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
-
+    def eval(cls, u):
         return ImmutableDenseMatrix([[dy(u[2]) - dz(u[1])],
                                      [dz(u[0]) - dx(u[2])],
                                      [dx(u[1]) - dy(u[0])]])
@@ -653,21 +604,15 @@ def eval(cls, *_args):
 #==============================================================================
 class Rot_2d(CalculusFunction):
 
-    nargs = None
+    nargs = None  # TODO [YG, 21.07.2022]: Check why we need this
     name = 'Grad'
 
-    def __new__(cls, *args, **options):
+    def __new__(cls, u, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
-            r = cls.eval(*args)
-        else:
-            r = None
-
-        if r is None:
-            return Basic.__new__(cls, *args, **options)
+            return cls.eval(u)
         else:
-            return r
+            return Basic.__new__(cls, u, **options)
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -677,33 +622,21 @@ def __getitem__(self, indices, **kw_args):
             return Indexed(self, indices, **kw_args)
 
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
-
-        return ImmutableDenseMatrix([[dy(u)],[-dx(u)]])
+    def eval(cls, u):
+        return ImmutableDenseMatrix([[dy(u)], [-dx(u)]])
 
 #==============================================================================
 class DivBasic(CalculusFunction):
 
-    nargs = None
+    nargs = None  # TODO [YG, 21.07.2022]: Check why we need this
     name = 'Div'
 
-    def __new__(cls, *args, **options):
+    def __new__(cls, u, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
-            r = cls.eval(*args)
-        else:
-            r = None
-
-        if r is None:
-            return Basic.__new__(cls, *args, **options)
+            return cls.eval(u)
         else:
-            return r
+            return Basic.__new__(cls, u, **options)
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -712,70 +645,78 @@ def __getitem__(self, indices, **kw_args):
         else:
             return Indexed(self, indices, **kw_args)
 
-class Div_1d(DivBasic):
-
-    @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
-        return dx(u[0])
+#    @classmethod
+#    def scalar_eval(cls, u):
+#        return sum(dxi(u[i]) for i, dxi in enumerate(cls._derivatives))
+#
+#    @classmethod
+#    def eval(cls, u):
+#        if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1:
+#            div = [cls.scalar_eval(u[:, j]) for j in range(u.shape[1])]
+#            return ImmutableDenseMatrix(div)
+#        else:
+#            return cls.scalar_eval(u)
 
-class Div_2d(DivBasic):
 
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
-        if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1:
-            div = []
-            for j in range(u.shape[1]):
-                div.append(dx(u[0,j])+dy(u[1,j]))
+    def eval(cls, u):
+        if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1:
+            div = [cls.eval(u[:, j]) for j in range(u.shape[1])]
             return ImmutableDenseMatrix(div)
+        else:
+            return sum(dxi(u[i]) for i, dxi in enumerate(cls._derivatives))
 
-        return dx(u[0]) + dy(u[1])
-
-class Div_3d(DivBasic):
 
-    @classmethod
-    def eval(cls, *_args):
+class Div_1d(DivBasic):
+    _derivatives = (dx, )
 
-        if not _args:
-            return
 
-        u = _args[0]
+class Div_2d(DivBasic):
+    _derivatives = (dx, dy)
 
-        if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1:
-            div = []
-            for j in range(u.shape[1]):
-                div.append(dx(u[0,j])+dy(u[1,j])+dz(u[2,j]))
-            return ImmutableDenseMatrix(div)
 
-        return dx(u[0]) + dy(u[1]) + dz(u[2])
+class Div_3d(DivBasic):
+    _derivatives = (dx, dy, dz)
+
+
+#class Div_1d(DivBasic):
+#    @classmethod
+#    def eval(cls, u):
+#        return dx(u[0])
+#
+#
+#class Div_2d(DivBasic):
+#    @classmethod
+#    def eval(cls, u):
+#        if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1:
+#            div = [dx(u[0, j]) + dy(u[1, j]) for j in range(u.shape[1])]
+#            return ImmutableDenseMatrix(div)
+#        else:
+#            return dx(u[0]) + dy(u[1])
+#
+#
+#class Div_3d(DivBasic):
+#    @classmethod
+#    def eval(cls, u):
+#        if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1:
+#            div = [dx(u[0, j]) + dy(u[1, j]) + dz(u[2, j])
+#                    for j in range(u.shape[1])]
+#            return ImmutableDenseMatrix(div)
+#        else:
+#            return dx(u[0]) + dy(u[1]) + dz(u[2])
 
 #==============================================================================
 class LaplaceBasic(CalculusFunction):
 
-    nargs = None
+    nargs = None  # TODO [YG, 21.07.2022]: Check why we need this
     name = 'Laplace'
 
-    def __new__(cls, *args, **options):
+    def __new__(cls, u, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
-            r = cls.eval(*args)
+            return cls.eval(u)
         else:
-            r = None
-
-        if r is None:
-            return Basic.__new__(cls, *args, **options)
-        else:
-            return r
+            return Basic.__new__(cls, u, **options)
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -784,66 +725,37 @@ def __getitem__(self, indices, **kw_args):
         else:
             return Indexed(self, indices, **kw_args)
 
-class Laplace_1d(LaplaceBasic):
-
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
+    def eval(cls, u):
         if isinstance(u, VectorFunction):
             raise NotImplementedError('TODO')
 
-        return dx(dx(u))
+        return sum(dxi(dxi(u)) for dxi in cls._derivatives)
 
-class Laplace_2d(LaplaceBasic):
 
-    @classmethod
-    def eval(cls, *_args):
+class Laplace_1d(LaplaceBasic):
+    _derivatives = (dx,)
 
-        if not _args:
-            return
 
-        u = _args[0]
-        if isinstance(u, VectorFunction):
-            raise NotImplementedError('TODO')
+class Laplace_2d(LaplaceBasic):
+    _derivatives = (dx, dy)
 
-        return dx(dx(u)) + dy(dy(u))
 
 class Laplace_3d(LaplaceBasic):
-
-    @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
-        if isinstance(u, VectorFunction):
-            raise NotImplementedError('TODO')
-
-        return dx(dx(u)) + dy(dy(u)) + dz(dz(u))
+    _derivatives = (dx, dy, dz)
 
 #==============================================================================
 class HessianBasic(CalculusFunction):
 
-    nargs = None
+    nargs = None  # TODO [YG, 21.07.2022]: Check why we need this
     name = 'Hessian'
 
-    def __new__(cls, *args, **options):
+    def __new__(cls, u, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
-            r = cls.eval(*args)
+            return cls.eval(u)
         else:
-            r = None
-
-        if r is None:
-            return Basic.__new__(cls, *args, **options)
-        else:
-            return r
+            return Basic.__new__(cls, u, **options)
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -852,50 +764,43 @@ def __getitem__(self, indices, **kw_args):
         else:
             return Indexed(self, indices, **kw_args)
 
-class Hessian_1d(HessianBasic):
+#    @classmethod
+#    def eval(cls, u):
+#        if isinstance(u, VectorFunction):
+#            raise NotImplementedError('TODO')
+#
+#        hessian = ImmutableDenseMatrix(
+#            [[di(dj(u)) for dj in cls._derivatives]
+#                        for di in cls._derivatives]
+#        )
+#
+#        return hessian[0, 0] if hessian.shape == (1, 1) else hessian
 
+class Hessian_1d(HessianBasic):
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
+    def eval(cls, u):
         if isinstance(u, VectorFunction):
             raise NotImplementedError('TODO')
+        return dx(u)
 
-        return dx(dx(u))
 
 class Hessian_2d(HessianBasic):
-
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
+    def eval(cls, u):
         if isinstance(u, VectorFunction):
             raise NotImplementedError('TODO')
-
         return ImmutableDenseMatrix([[dx(dx(u)), dx(dy(u))],
-                       [dx(dy(u)), dy(dy(u))]])
+                                     [dx(dy(u)), dy(dy(u))]])
 
-class Hessian_3d(HessianBasic):
 
+class Hessian_3d(HessianBasic):
     @classmethod
-    def eval(cls, *_args):
-
-        if not _args:
-            return
-
-        u = _args[0]
+    def eval(cls, u):
         if isinstance(u, VectorFunction):
             raise NotImplementedError('TODO')
-
         return ImmutableDenseMatrix([[dx(dx(u)), dx(dy(u)), dx(dz(u))],
-                       [dx(dy(u)), dy(dy(u)), dy(dz(u))],
-                       [dx(dz(u)), dy(dz(u)), dz(dz(u))]])
+                                     [dx(dy(u)), dy(dy(u)), dy(dz(u))],
+                                     [dx(dz(u)), dy(dz(u)), dz(dz(u))]])
 
 #==============================================================================
 class BracketBasic(CalculusFunction):