diff --git a/hysop/backend/device/opencl/operator/derivative.py b/hysop/backend/device/opencl/operator/derivative.py
index f7424484377cb9acbaf11f70a77c23d34626f21e..c074092fb71580ba372a2510f4500784863d0714 100644
--- a/hysop/backend/device/opencl/operator/derivative.py
+++ b/hysop/backend/device/opencl/operator/derivative.py
@@ -19,7 +19,8 @@ from hysop.operator.base.custom_symbolic_operator import SymbolicExpressionParse
 from hysop.symbolic.relational import Assignment
 
 
-class OpenClFiniteDifferencesSpaceDerivative(FiniteDifferencesSpaceDerivativeBase, OpenClSymbolic):
+class OpenClFiniteDifferencesSpaceDerivative(FiniteDifferencesSpaceDerivativeBase, 
+                                             OpenClSymbolic):
 
     @debug
     def __init__(self, **kwds):
@@ -71,6 +72,18 @@ class OpenClSpectralSpaceDerivative(SpectralSpaceDerivativeBase, OpenClSymbolic)
             Base class arguments.
         """
         super(OpenClSpectralSpaceDerivative, self).__init__(**kwds)
+        Fs  = self.Fin.s()
+        dFs = self.Fout.s()
+        if self.scale_by_field:
+            As = self.A.s()
+        else:
+            As = None
+        print
+        # TODO
+        for wn in self.tg._wave_numbers:
+            print wn
+            import sys 
+            sys.exit(1)
     
     @debug
     def setup(self, work):
@@ -92,7 +105,6 @@ class OpenClSpectralSpaceDerivative(SpectralSpaceDerivativeBase, OpenClSymbolic)
         self.dFout.exchange_ghosts()
 
     def compute_derivative(self):
-        from hysop.constants import BoxBoundaryCondition
         for nd_dkd in self.nd_dkds:
             self.Ft.output_buffer[...] *= nd_dkd
 
@@ -108,3 +120,4 @@ class OpenClSpectralSpaceDerivative(SpectralSpaceDerivativeBase, OpenClSymbolic)
                 out[...] *= scale()
         elif self.scale_by_value:
             out[...] *= scale
+
diff --git a/hysop/backend/device/opencl/operator/min_max.py b/hysop/backend/device/opencl/operator/min_max.py
index 0211260a69b48200b1bc68ad93975236dd1adee5..03a88e9e53e1c0eed3055a416b218df0dcd75292 100644
--- a/hysop/backend/device/opencl/operator/min_max.py
+++ b/hysop/backend/device/opencl/operator/min_max.py
@@ -4,7 +4,8 @@ from hysop.core.graph.graph import op_apply
 from hysop.operator.base.min_max import MinMaxFieldStatisticsBase, \
                                         MinMaxDerivativeStatisticsBase
 from hysop.backend.device.opencl.opencl_operator import OpenClOperator
-from hysop.backend.device.opencl.operator.derivative import OpenClSpaceDerivative
+from hysop.backend.device.opencl.operator.derivative import OpenClSpectralSpaceDerivative, \
+                                                            OpenClFiniteDifferencesSpaceDerivative
 
 
 class OpenClMinMaxFieldStatistics(MinMaxFieldStatisticsBase, OpenClOperator):
@@ -17,11 +18,21 @@ class OpenClMinMaxFieldStatistics(MinMaxFieldStatisticsBase, OpenClOperator):
         self.compute_statistics(**kwds)
 
 
-class OpenClMinMaxDerivativeStatistics(MinMaxDerivativeStatisticsBase, 
-                                               OpenClSpaceDerivative):
-    """OpenCl implementation backend of operator MinMaxDerivativeStatistics."""
+class OpenClMinMaxSpectralDerivativeStatistics(MinMaxDerivativeStatisticsBase, 
+                                               OpenClSpectralSpaceDerivative):
+    """OpenCl implementation backend of operator MinMaxSpectralDerivativeStatistics."""
     @op_apply
     def apply(self, **kwds):
         """Compute derivative and than statistics."""
-        super(OpenClMinMaxDerivativeStatistics, self).apply(**kwds)
+        super(OpenClMinMaxSpectralDerivativeStatistics, self).apply(**kwds)
+        self.compute_statistics(**kwds)
+
+
+class OpenClMinMaxFiniteDifferencesDerivativeStatistics(MinMaxDerivativeStatisticsBase, 
+                                                        OpenClFiniteDifferencesSpaceDerivative):
+    """OpenCl implementation backend of operator MinMaxFiniteDifferencesDerivativeStatistics."""
+    @op_apply
+    def apply(self, **kwds):
+        """Compute derivative and than statistics."""
+        super(OpenClMinMaxFiniteDifferencesDerivativeStatistics, self).apply(**kwds)
         self.compute_statistics(**kwds)
diff --git a/hysop/backend/host/python/operator/derivative.py b/hysop/backend/host/python/operator/derivative.py
index c265eb8c5e65ad63b6c6f4be2f38659fe74cc0ab..48c593461034d14246a0f7de000a8f4decb20b71 100644
--- a/hysop/backend/host/python/operator/derivative.py
+++ b/hysop/backend/host/python/operator/derivative.py
@@ -5,9 +5,8 @@ from hysop.operator.base.derivative import FiniteDifferencesSpaceDerivativeBase,
 from hysop.backend.host.host_operator import HostOperator
 from hysop.tools.decorators import debug
 from hysop.core.graph.graph import op_apply
-from hysop.numerics.stencil.stencil_generator import StencilGenerator, CenteredStencilGenerator, MPQ
-
-
+from hysop.numerics.stencil.stencil_generator import StencilGenerator, \
+        CenteredStencilGenerator, MPQ
 
 class PythonSpectralSpaceDerivative(SpectralSpaceDerivativeBase, HostOperator):
     """
@@ -107,7 +106,8 @@ class PythonFiniteDifferencesSpaceDerivative(FiniteDifferencesSpaceDerivativeBas
         d = self.d
         
         # set min_ghosts for input field
-        requirements = super(PythonFiniteDifferencesSpaceDerivative, self).get_field_requirements()
+        requirements = super(PythonFiniteDifferencesSpaceDerivative, 
+                self).get_field_requirements()
         for is_input, (field, td, req) in requirements.iter_requirements():
             if (field is self.Fin):
                 ghosts = req.min_ghosts.copy()
diff --git a/hysop/backend/host/python/operator/min_max.py b/hysop/backend/host/python/operator/min_max.py
index f07b16951a182eac826d1900c26da6bf9dc162e2..0f7ce42cad1887b1fabd0e0cc68bcee6a55ec1c9 100644
--- a/hysop/backend/host/python/operator/min_max.py
+++ b/hysop/backend/host/python/operator/min_max.py
@@ -3,9 +3,11 @@ from hysop.core.graph.graph import op_apply
 from hysop.operator.base.min_max import MinMaxFieldStatisticsBase, \
                                         MinMaxDerivativeStatisticsBase
 from hysop.backend.host.host_operator import HostOperator
-from hysop.backend.host.python.operator.derivative import PythonSpaceDerivative
+from hysop.backend.host.python.operator.derivative import PythonSpectralSpaceDerivative, \
+        PythonFiniteDifferencesSpaceDerivative
 
-class PythonMinMaxFieldStatistics(MinMaxFieldStatisticsBase, HostOperator):
+class PythonMinMaxFieldStatistics(MinMaxFieldStatisticsBase, 
+                                  HostOperator):
     """Python implementation backend of operator MinMaxFieldStatistics."""
 
     @debug
@@ -20,11 +22,21 @@ class PythonMinMaxFieldStatistics(MinMaxFieldStatisticsBase, HostOperator):
         self.compute_statistics(**kwds)
 
 
-class PythonMinMaxDerivativeStatistics(MinMaxDerivativeStatisticsBase, 
-                                               PythonSpaceDerivative):
-    """Python implementation backend of operator MinMaxDerivativeStatistics."""
+class PythonMinMaxSpectralDerivativeStatistics(MinMaxDerivativeStatisticsBase, 
+                                               PythonSpectralSpaceDerivative):
+    """Python implementation backend of operator MinMaxSpectralDerivativeStatistics."""
     @op_apply
     def apply(self, **kwds):
-        """Compute derivative and than statistics."""
-        super(PythonMinMaxDerivativeStatistics, self).apply(**kwds)
+        """Compute derivative and then statistics."""
+        super(PythonMinMaxSpectralDerivativeStatistics, self).apply(**kwds)
+        self.compute_statistics(**kwds)
+
+
+class PythonMinMaxFiniteDifferencesDerivativeStatistics(MinMaxDerivativeStatisticsBase, 
+                                                        PythonFiniteDifferencesSpaceDerivative):
+    """Python implementation backend of operator MinMaxFiniteDifferencesDerivativeStatistics."""
+    @op_apply
+    def apply(self, **kwds):
+        """Compute derivative and then statistics."""
+        super(PythonMinMaxFiniteDifferencesDerivativeStatistics, self).apply(**kwds)
         self.compute_statistics(**kwds)
diff --git a/hysop/fields/default_fields.py b/hysop/fields/default_fields.py
index 2710608f8ab8de65499495957861d3eabe92228c..39d9f2db563394553648a2fe66a1ba870211c02f 100644
--- a/hysop/fields/default_fields.py
+++ b/hysop/fields/default_fields.py
@@ -1,36 +1,61 @@
-from hysop.tools.types import first_not_None
+from hysop.tools.types import first_not_None, check_instance
 from hysop.tools.sympy_utils import greak, Greak, subscripts
-from hysop.fields.continuous_field import Field
+from hysop.fields.continuous_field import Field, TensorField 
+from hysop.tools.numpywrappers import npw
+from hysop.constants import BoxBoundaryCondition, BoundaryCondition
 
-def VelocityField(domain, name=None, pretty_name=None, 
-                    is_vector=True, **kwds):
+
+def VelocityField(domain, name=None, pretty_name=None, **kwds):
     name        = first_not_None(name, 'U')
     pretty_name = first_not_None(pretty_name, greak[20])
-    is_vector   = first_not_None(is_vector, True)
-    return Field(domain=domain, name=name, pretty_name=pretty_name, 
-            is_vector=is_vector, **kwds)
+    lboundaries, rboundaries = domain.lboundaries, domain.rboundaries
+    dim = domain.dim
+    def velocity_boundaries(boundaries, component):
+        check_instance(boundaries, npw.ndarray, values=BoxBoundaryCondition)
+        fboundaries = npw.empty_like(boundaries)
+        fboundaries[...] = None
+        for (i,bd) in enumerate(boundaries):
+            is_normal = (dim-i-1==component)
+            if (bd is BoxBoundaryCondition.PERIODIC):
+                fbd = BoundaryCondition.PERIODIC
+            elif (bd is BoxBoundaryCondition.SYMMETRIC): # (normal to boundary velocity = 0)
+                if is_normal:
+                    fbd = BoundaryCondition.HOMOGENEOUS_DIRICHLET
+                else:
+                    fbd = BoundaryCondition.HOMOGENEOUS_NEUMANN
+            elif (bd is BoxBoundaryCondition.OUTFLOW): # (velocity is normal to boundary)
+                if is_normal:
+                    fbd = BoundaryCondition.HOMOGENEOUS_NEUMANN
+                else:
+                    fbd = BoundaryCondition.HOMOGENEOUS_DIRICHLET
+            else:
+                msg='FATAL ERROR: Unknown domain boundary condition {}.'
+                msg=msg.format(bd)
+                raise NotImplementedError(msg)
+            fboundaries[i] = fbd
+        check_instance(fboundaries, npw.ndarray, values=BoundaryCondition)
+        return fboundaries
+    def _make_field(idx, **fkwds): 
+        # Adapt velocity boundaries to domain boundaries
+        component, = idx
+        fkwds['lboundaries'] = velocity_boundaries(lboundaries, component)
+        fkwds['rboundaries'] = velocity_boundaries(rboundaries, component)
+        return TensorField.default_make_field(idx=idx, **fkwds)
+    kwds.setdefault('make_field', _make_field)
+    kwds.setdefault('is_vector', True)
+    return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
 
-def VorticityField(domain, name=None, pretty_name=None, 
-                    nb_components=None, **kwds):
+
+def VorticityField(velocity, name=None, pretty_name=None, **kwds):
+    # vorticity domain domain boundaries are deduced from velocity boundary conditions
+    check_instance(velocity, Field)
+    domain = velocity.domain
+    assert velocity.nb_components == domain.dim, 'Invalid velocity Field.'
     name        = first_not_None(name, 'W')
     pretty_name = first_not_None(pretty_name, greak[24])
-    if (nb_components is None):
-        if (domain.dim == 2):
-            nb_components = 1
-        elif (domain.dim == 3):
-            nb_components = 3
-        else:
-            msg='Cannot deduce the number of components of the vorticity '
-            msg+='for a {}d domain.'.format(domain.dim)
-            raise ValueError(msg)
-    return Field(domain=domain, name=name, pretty_name=pretty_name, 
-                    nb_components=nb_components, **kwds)
-
-    
-def CurvatureField(domain, name=None, pretty_name=None, **kwds):
-    name        = first_not_None(name, 'kappa')
-    pretty_name = first_not_None(pretty_name, greak[9])
-    return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
+    return velocity.curl(name=name, pretty_name=pretty_name, 
+            scalar_name_prefix=name, scalar_pretty_name_prefix=pretty_name,
+            **kwds)
 
 
 def DensityField(domain, name=None, pretty_name=None, **kwds):
@@ -48,12 +73,21 @@ def ViscosityField(domain, name=None, pretty_name=None, mu=False, **kwds):
         pretty_name = first_not_None(pretty_name, greak[12])
     return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
 
+
 def LevelSetField(domain, name=None, pretty_name=None, **kwds):
     name = first_not_None(name, 'phi')
     pretty_name = first_not_None(pretty_name, Greak[21])
     return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
 
+
 def PenalizationField(domain, name=None, pretty_name=None, **kwds):
     name = first_not_None(name, 'lambda')
     pretty_name = first_not_None(pretty_name, greak[10])
     return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
+
+
+def CurvatureField(domain, name=None, pretty_name=None, **kwds):
+    name        = first_not_None(name, 'kappa')
+    pretty_name = first_not_None(pretty_name, greak[9])
+    return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
+
diff --git a/hysop/operator/base/derivative.py b/hysop/operator/base/derivative.py
index d838dccf20d1e4dbc64c145a259f24bfb9f4a235..ee80ce6c98f713bd881de3136e02074d28a5527d 100644
--- a/hysop/operator/base/derivative.py
+++ b/hysop/operator/base/derivative.py
@@ -249,7 +249,6 @@ class SpectralSpaceDerivativeBase(SpectralOperatorBase, SpaceDerivativeBase):
         Bt = tg.require_backward_transform(dF, axes=axes)
         assert (Ft.output_dtype == Bt.input_dtype)
         
-
         dFt = Ft.s
         assert len(derivative)==F.dim
         for (i,di) in enumerate(derivative):
diff --git a/hysop/operator/base/spectral_operator.py b/hysop/operator/base/spectral_operator.py
index d616b63b2876af8427d8b5c7506d940de2c91e79..b0990f8e3d318230693109ef8b244102d9371673 100644
--- a/hysop/operator/base/spectral_operator.py
+++ b/hysop/operator/base/spectral_operator.py
@@ -474,7 +474,7 @@ class SpectralTransformGroup(object):
             *DCT-III -> DCT-II
             *DST-I   -> DST-I
             *DST-III -> DST-II
-        else, if order[axis] is odd:
+        else: (if order[axis] is odd)
             *C2C     -> C2C
             *R2C     -> C2R
             *DCT-I   -> DST-I
diff --git a/hysop/operator/derivative.py b/hysop/operator/derivative.py
index 54de8b698b998721763b5339c2ca8d24ef213b6e..b57a7a6ca7adad5d6b6e1524333028bc015ac241 100644
--- a/hysop/operator/derivative.py
+++ b/hysop/operator/derivative.py
@@ -193,8 +193,10 @@ class MultiSpaceDerivatives(DirectionalOperatorGeneratorI, ComputationalGraphNod
 
          Refer to hysop.operator.derivative.SpaceDerivative for all arguments.
          """
-         if not issubclass(cls, SpaceDerivative) or (cls is SpaceDerivative):
-             msg="cls should be a subclass of SpaceDerivative, got {}."
+         from hysop.operator.min_max import MinMaxDerivativeStatistics
+         if not issubclass(cls, (SpaceDerivative, MinMaxDerivativeStatistics)) or \
+                (cls in (SpaceDerivative, MinMaxDerivativeStatistics)):
+             msg="cls should be a subclass of SpaceDerivative or MinMaxSpaceDerivativeStatistics, got {}."
              msg+='\ncls MRO is:\n  '
              msg+='\n  '.join(str(t) for t in cls.__mro__)
              msg=msg.format(cls)
@@ -384,9 +386,11 @@ class Gradient(MultiSpaceDerivatives):
         base_kwds.update(dict(
                 candidate_input_tensors=(F,),
                 candidate_output_tensors=(gradF,)))
-
-        if not issubclass(cls, SpaceDerivative) or (cls is SpaceDerivative):
-            msg="cls should be a subclass of SpaceDerivative, got {}."
+        
+        from hysop.operator.min_max import MinMaxDerivativeStatistics
+        if not issubclass(cls, (SpaceDerivative, MinMaxDerivativeStatistics)) or \
+                (cls in (SpaceDerivative, MinMaxDerivativeStatistics)):
+            msg="cls should be a subclass of SpaceDerivative or MinMaxSpaceDerivativeStatistics, got {}."
             msg+='\ncls MRO is:\n  '
             msg+='\n  '.join(str(t) for t in cls.__mro__)
             msg=msg.format(cls)
diff --git a/hysop/operator/min_max.py b/hysop/operator/min_max.py
index 38814c859abe86e19805e436ee68940d3f6cd2ab..57ffae91c5d348666b10440b711810abe1230bed 100644
--- a/hysop/operator/min_max.py
+++ b/hysop/operator/min_max.py
@@ -1,8 +1,10 @@
 """
 @file min_max.py
-MinMaxFieldStatistics:      compute min(f), max(f) and/or max(|f|) for a given field f, component-wise. 
-MinMaxDerivativeStatistics: compute min(d^k(Fi)/dXj^k), max(d^kFi/dXj^k) and/or max(|dFi/dXj|) for a given field, component, direction and order.
-MinMaxGradientStatistics:   compute min(dFi/dXj), max(dFi/dXj) and/or max(|dFi/dXj|) for a given field, up to all components in all directions.
+MinMaxFieldStatistics: compute min(f), max(f) and/or max(|f|) for a given field f.
+MinMaxDerivativeStatistics: compute min(d^k(Fi)/dXj^k), max(d^kFi/dXj^k) and/or max(|dFi/dXj|)
+                            for a given field, component, direction and order.
+MinMaxGradientStatistics: compute min(dFi/dXj), max(dFi/dXj) and/or max(|dFi/dXj|) 
+                          for a given field, up to all components in all directions.
 """
 from hysop import vprint
 from hysop.constants import Backend, Implementation
@@ -53,7 +55,8 @@ class MinMaxFieldStatistics(ComputationalGraphNodeFrontend):
         MinMaxFieldStatistics can compute some commonly required Field statistics:
             Fmin:  component-wise min values of the field.
             Fmax:  component-wise max values of the field.
-            Finf:  component-wise max values of the absolute value of the field (computed using Fmin and Fmax).
+            Finf:  component-wise max values of the absolute value of the field (computed using
+                                                                                 Fmin and Fmax).
 
         All statistics are only computed if explicitely requested by user,
           unless required to compute another user-required statistic, see Notes.
@@ -197,8 +200,8 @@ class MinMaxDerivativeStatistics(ComputationalGraphNodeFrontend):
                     derivative of the field (computed using Fmin and Fmax).
         
         First compute the derivative of a component of a field F in a given direction
-        at a given order and on a given backend out of place in a specific output component of dF. 
-        The derivative is then possibly scaled by another field/parameter/value A.
+        at a given order and on a given backend out of place in a specific output component of 
+        dF. The derivative is then possibly scaled by another field/parameter/value A.
 
         After the scaled derivative has been computed, compute user requested statistics
         (min and max values) on this new field and scale those statistics by other scaling 
@@ -318,7 +321,8 @@ class MinMaxDerivativeStatistics(ComputationalGraphNodeFrontend):
         check_instance(direction, int, allow_none=True)
         check_instance(out_component, int, allow_none=True)
         check_instance(coeffs, dict, keys=str, values=(int, float, npw.number), allow_none=True)
-        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors, allow_none=True)
+        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors,
+                        allow_none=True)
         check_instance(name, str, allow_none=True)
         check_instance(pbasename, str, allow_none=True)
         check_instance(ppbasename, (str, unicode), allow_none=True)
@@ -353,10 +357,16 @@ class MinMaxDerivativeStatistics(ComputationalGraphNodeFrontend):
 
 
 class MinMaxSpectralDerivativeStatistics(MinMaxDerivativeStatistics):
+    """
+    Operator frontend to compute min and max statistics on a specific 
+    derivative of a field component using the spectral method.
+    """
     @classmethod
     def implementations(cls):
-        from hysop.backend.host.python.operator.min_max   import PythonMinMaxSpectralDerivativeStatistics
-        from hysop.backend.device.opencl.operator.min_max import OpenClMinMaxSpectralDerivativeStatistics
+        from hysop.backend.host.python.operator.min_max import \
+                PythonMinMaxSpectralDerivativeStatistics
+        from hysop.backend.device.opencl.operator.min_max import \
+                OpenClMinMaxSpectralDerivativeStatistics
         implementations = {
                 Implementation.PYTHON: PythonMinMaxSpectralDerivativeStatistics,
                 Implementation.OPENCL: OpenClMinMaxSpectralDerivativeStatistics
@@ -369,10 +379,16 @@ class MinMaxSpectralDerivativeStatistics(MinMaxDerivativeStatistics):
 
 
 class MinMaxFiniteDifferencesDerivativeStatistics(MinMaxDerivativeStatistics):
+    """
+    Operator frontend to compute min and max statistics on a specific 
+    derivative of a field component using finite differences.
+    """
     @classmethod
     def implementations(cls):
-        from hysop.backend.host.python.operator.min_max   import PythonMinMaxFiniteDifferencesDerivativeStatistics
-        from hysop.backend.device.opencl.operator.min_max import OpenClMinMaxFiniteDifferencesDerivativeStatistics
+        from hysop.backend.host.python.operator.min_max import \
+                PythonMinMaxFiniteDifferencesDerivativeStatistics
+        from hysop.backend.device.opencl.operator.min_max import \
+                OpenClMinMaxFiniteDifferencesDerivativeStatistics
         implementations = {
                 Implementation.PYTHON: PythonMinMaxFiniteDifferencesDerivativeStatistics,
                 Implementation.OPENCL: OpenClMinMaxFiniteDifferencesDerivativeStatistics
@@ -397,6 +413,7 @@ class MinMaxGradientStatistics(Gradient):
             all_quiet=True, print_tensors=True,
             name=None, pretty_name=None, pbasename=None, ppbasename=None,
             variables=None, implementation=None, base_kwds=None, 
+            cls=MinMaxFiniteDifferencesDerivativeStatistics,
             **kwds):
         """
         Create an operator generator that yields a sequence of operators
@@ -408,38 +425,30 @@ class MinMaxGradientStatistics(Gradient):
             Finf:  component-wise and direction-wise max values of the absolute value of the 
                         gradient of the field (computed using Fmin and Fmax).
 
-        Derivatives can be computed with respect to specific directions and not necessarily in all directions.
-        To restrict the number of components, take a tensor view on F (and gradF).
-        
-        ==============================================================================================================
-                   dF0/dX0 ... dF0/dXn                             0    ...    n
-                      .     .     .                                .     .     .
-        grad(F) =     .     .     .     are mapped to components   .     .     .          in the output statistics.
-                      .     .     .                                .     .     .             n = nb_directions - 1
-                   dFm/dX0 ... dFm/dXn                          m(n+1)   .  (m+1)(n+1)-1     m = F.nb_components - 1
-        ==============================================================================================================
+        Derivatives can be computed with respect to specific directions and not necessarily in 
+        all directions. To restrict the number of components, take a tensor view on F (and gradF).
         
+        ----------------------------------------------
         Let k = idx + (j,)
-            --------------------
-            gradF[k] = dF[idx]/dXd
-            --------------------
-            Fmax[k]  = Smin * min( dF[idx]/dXd )
-            Fmin[k]  = Smax * max( dF[idx]/dXd )
-            Finf[k]  = Sinf * max( |Fmin[k]|, |Fmax[k]| )
-
-            where F is an input field,
-                  nb_components = F.nb_components = np.prod(F.shape)
-                  nb_directions = min( F.dim, len(directions))
-                  gradF is an optional output tensor field,
-                  idx is contained in numpy.ndindex(*F.shape)
-                  0 <= j < nb_directions
-                  d = directions[j]
-                  Fmin = created or supplied TensorParameter of shape F.shape + (nb_directions,).
-                  Fmax = created or supplied TensorParameter of shape F.shape + (nb_directions,).
-                  Finf = created or supplied TensorParameter of shape F.shape + (nb_directions,).
-                  Smin = coeffs['Fmin'] or coeffs['Fmin'][k]
-                  Smax = coeffs['Fmax'] or coeffs['Fmax'][k]
-                  Sinf = coeffs['Finf'] or coeffs['Fmax'][k]
+        gradF[k] = dF[idx]/dXd
+        ----------------------------------------------
+        Fmax[k]  = Smin * min( dF[idx]/dXd )
+        Fmin[k]  = Smax * max( dF[idx]/dXd )
+        Finf[k]  = Sinf * max( |Fmin[k]|, |Fmax[k]| )
+        ----------------------------------------------
+        where F is an input field,
+              nb_components = F.nb_components = np.prod(F.shape)
+              nb_directions = min( F.dim, len(directions))
+              gradF is an optional output tensor field,
+              idx is contained in numpy.ndindex(*F.shape)
+              0 <= j < nb_directions
+              d = directions[j]
+              Fmin = created or supplied TensorParameter of shape F.shape + (nb_directions,).
+              Fmax = created or supplied TensorParameter of shape F.shape + (nb_directions,).
+              Finf = created or supplied TensorParameter of shape F.shape + (nb_directions,).
+              Smin = coeffs['Fmin'] or coeffs['Fmin'][k]
+              Smax = coeffs['Fmax'] or coeffs['Fmax'][k]
+              Sinf = coeffs['Finf'] or coeffs['Fmax'][k]
 
         All statistics are only computed if explicitely required by user,
           unless required to compute another required statistic, see Notes.
@@ -521,7 +530,8 @@ class MinMaxGradientStatistics(Gradient):
         check_instance(directions, tuple, values=int, allow_none=True, 
                 minval=0, maxval=F.dim-1, minsize=1, unique=True)
         check_instance(coeffs, dict, keys=str, values=(int, float, npw.number), allow_none=True)
-        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors, allow_none=True)
+        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors, 
+                        allow_none=True)
         check_instance(name, str, allow_none=True)
         check_instance(pbasename, str, allow_none=True)
         check_instance(ppbasename, (str,unicode), allow_none=True)
@@ -572,7 +582,8 @@ class MinMaxGradientStatistics(Gradient):
         ppbasename = first_not_None(ppbasename, gradF.pretty_name)
         
         names = { k: v.format(pbasename) for (k,v) in _names.iteritems() }
-        pretty_names = { k: v.format(ppbasename.decode('utf-8')) for (k,v) in _pretty_names.iteritems() }
+        pretty_names = { k: v.format(ppbasename.decode('utf-8')) 
+                            for (k,v) in _pretty_names.iteritems() }
 
         def make_param(k, quiet):
             return TensorParameter(name=names[k], pretty_name=pretty_names[k],
@@ -621,11 +632,12 @@ class MinMaxGradientStatistics(Gradient):
                 stats = extra_params.setdefault(statname, gradF.new_empty_array())
                 stats[idx] = S
             extra_params['name'][idx] = name.format(Fi.name)
-            extra_params['pretty_name'][idx] = pretty_name.format(Fi.pretty_name.decode('utf-8')).encode('utf-8')
+            extra_params['pretty_name'][idx] = pretty_name.format(
+                                                Fi.pretty_name.decode('utf-8')).encode('utf-8')
 
-        super(MinMaxGradientStatistics, self).__init__(F=F, gradF=gradF, 
+        super(MinMaxGradientStatistics, self).__init__(F=F, gradF=gradF,
                 directions=directions, extra_params=extra_params,
-                variables=variables, **kwds)
+                cls=cls, variables=variables, **kwds)
 
         # add a phony operator to gather parameter views
         class MergeTensorViewsOperator(ComputationalGraphOperator):
diff --git a/hysop/operators.py b/hysop/operators.py
index 86aa1220ef7be0f26e338b49a50701eb69579a63..f113d34531df35662f1b4e59b676d701c9e4a0b7 100644
--- a/hysop/operators.py
+++ b/hysop/operators.py
@@ -28,7 +28,8 @@ from hysop.operator.derivative import SpaceDerivative,                  \
                                       Gradient
 
 from hysop.operator.min_max import MinMaxFieldStatistics,      \
-                                   MinMaxDerivativeStatistics, \
+                                   MinMaxFiniteDifferencesDerivativeStatistics, \
+                                   MinMaxSpectralDerivativeStatistics, \
                                    MinMaxGradientStatistics
 
 from hysop.numerics.splitting.strang import StrangSplitting
diff --git a/requirements.txt b/requirements.txt
index 502b8de768a264fa733d349d0d466826316bdf42..6c38fcf73cdaee2291bee0527459be7cdaac4555 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -15,6 +15,5 @@ argparse_color_formatter
 primefac
 pyopencl
 pyfftw
-gpyfft
 mpi4py
 matplotlib