From 7ff04dc9bd7b96806ba8791d261758c5d6fa54ae Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Wed, 22 May 2019 18:44:35 +0200
Subject: [PATCH] fixed everything for ci

---
 examples/analytic/analytic.py                 | 20 ++--
 examples/scalar_advection/levelset.py         | 21 ++--
 examples/taylor_green/taylor_green.py         |  2 +-
 .../device/opencl/operator/analytic.py        |  6 +-
 .../backend/host/python/operator/analytic.py  | 15 +--
 .../operator/directional/advection_dir.py     |  7 +-
 hysop/core/graph/computational_node.py        |  6 ++
 hysop/fields/cartesian_discrete_field.py      | 17 ++--
 hysop/fields/tests/test_cartesian.py          | 22 ++---
 hysop/methods.py                              |  2 +-
 hysop/numerics/tests/test_fft.py              |  5 +-
 hysop/operator/analytic.py                    | 99 ++++++++++++++++++-
 hysop/operator/base/spatial_filtering.py      | 13 ++-
 hysop/operator/spatial_filtering.py           | 17 ++--
 hysop/operator/tests/test_absorption.py       | 13 +--
 hysop/operator/tests/test_analytic.py         | 26 +++--
 hysop/operator/tests/test_diffusion.py        | 20 ++--
 .../tests/test_directional_advection.py       |  2 +-
 .../tests/test_directional_diffusion.py       | 19 ++--
 .../tests/test_directional_stretching.py      | 12 +--
 hysop/operator/tests/test_fd_derivative.py    | 13 +--
 hysop/operator/tests/test_poisson.py          | 22 ++---
 hysop/operator/tests/test_poisson_curl.py     | 23 ++---
 .../operator/tests/test_restriction_filter.py | 19 ++--
 .../tests/test_solenoidal_projection.py       | 28 +++---
 hysop/operator/tests/test_spectral_curl.py    | 26 ++---
 .../tests/test_spectral_derivative.py         | 22 ++---
 hysop/operator/tests/test_transpose.py        | 21 ++--
 hysop/testsenv.py                             | 18 ++--
 hysop/topology/cartesian_descriptor.py        | 18 +++-
 30 files changed, 315 insertions(+), 239 deletions(-)

diff --git a/examples/analytic/analytic.py b/examples/analytic/analytic.py
index 09f01a5cd..569583479 100755
--- a/examples/analytic/analytic.py
+++ b/examples/analytic/analytic.py
@@ -6,10 +6,10 @@ def compute(args):
     '''
     HySoP Analytic Example: Initialize a field with a space and time dependent analytic formula.
     '''
-    from hysop import Field, Box, MPIParams, \
+    from hysop import Field, Box, IOParams, MPIParams, \
                       Simulation, Problem, ScalarParameter
     from hysop.constants import Implementation
-    from hysop.operators import AnalyticField
+    from hysop.operators import AnalyticField, HDF_Writer
     
     # Define domain
     npts = args.npts
@@ -57,10 +57,10 @@ def compute(args):
         # tuple 'data'. Coordinates will be passed as a tuple as a second 
         # argument. Finally extra arguments (here t) are passed last.
         # Note that t is a ScalarParameter, so we evaluate it to get its value.
-        def compute_scalar(data, coords, t):
-            data[0][...] = (1.0/(1.0+0.1*t()))
-            for x in coords[0]:
-                data[0][...] *= np.cos(x-t())
+        def compute_scalar(data, coords, component, t):
+            data[...] = (1.0/(1.0+0.1*t()))
+            for x in coords:
+                data[...] *= np.cos(x-t())
     elif (impl is Implementation.OPENCL):
         # With the opencl codegen implementation we use a symbolic expression
         # generated using sympy. OpenCL code will be automatically generated,
@@ -83,11 +83,15 @@ def compute(args):
                 **op_kwds)
 
     # Write output field at given frequency
-    analytic.dump_outputs(fields=scalar, frequency=args.dump_freq, filename='F', **op_kwds)
+    io_params = IOParams(filename='analytic', frequency=args.dump_freq)
+    df = HDF_Writer(name='S',
+                     io_params=io_params,
+                     variables={scalar: npts},
+                     **op_kwds)
    
     # Create the problem we want to solve and insert our operator
     problem = Problem()
-    problem.insert(analytic)
+    problem.insert(analytic, df)
     problem.build(args)
 
     # If a visu_rank was provided, and show_graph was set,
diff --git a/examples/scalar_advection/levelset.py b/examples/scalar_advection/levelset.py
index 287974551..1681b1b13 100644
--- a/examples/scalar_advection/levelset.py
+++ b/examples/scalar_advection/levelset.py
@@ -91,17 +91,16 @@ def compute(args):
     # Setup implementation specific variables
     if (impl is Implementation.PYTHON) or (impl is Implementation.FORTRAN):
         sin, cos = np.sin, np.cos
-        def compute_velocity(data, coords, t):
-            (x,y,z) = coords[0]
-            data[0][...] = 2.*sin(pi*x)**2*sin(2*pi*y)*sin(2.*pi*z)*cos(t()*pi/3.)
-            (x,y,z) = coords[1]
-            data[1][...] = -sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)*cos(t()*pi/3.)
-            (x,y,z) = coords[2]
-            data[2][...] = -sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2*cos(t()*pi/3.)
-        def compute_norm(data, coords):
-            data[0][...] = np.sqrt(data[0]**2 + data[1]**2 + data[2]**2)
-        def compute_volume(data, coords, S):
-            data[0][...] = (S[0] <= 0.5)
+        def compute_velocity(data, coords, component, t):
+            (x,y,z) = coords
+            fn = [ +2.*sin(pi*x)**2*sin(2*pi*y)*sin(2.*pi*z)*cos(t()*pi/3.),
+                   -sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)*cos(t()*pi/3.),
+                   -sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2*cos(t()*pi/3.)]
+            data[...] = fn[component]
+        def compute_norm(data, coords, component):
+            data[...] = np.sqrt(data[0]**2 + data[1]**2 + data[2]**2)
+        def compute_volume(data, coords, component, S):
+            data[...] = (S[0] <= 0.5)
     elif (impl is Implementation.OPENCL):
         from hysop.symbolic.relational import LogicalLE
         sin, cos = sm.sin, sm.cos
diff --git a/examples/taylor_green/taylor_green.py b/examples/taylor_green/taylor_green.py
index b470807c0..c00e88b98 100644
--- a/examples/taylor_green/taylor_green.py
+++ b/examples/taylor_green/taylor_green.py
@@ -10,7 +10,7 @@ cos = np.cos
 sin = np.sin
 
 ## Function to compute initial vorticity
-def init_vorticity(data, coords, component=None):
+def init_vorticity(data, coords, component):
     # Ux = sin(x) * cos(y) * cos(z)
     # Uy = - cos(x) * sin(y) * cos(z)
     # Uz = 0
diff --git a/hysop/backend/device/opencl/operator/analytic.py b/hysop/backend/device/opencl/operator/analytic.py
index d490f3654..916ea10ae 100644
--- a/hysop/backend/device/opencl/operator/analytic.py
+++ b/hysop/backend/device/opencl/operator/analytic.py
@@ -2,7 +2,7 @@
 from hysop.deps import sm
 from hysop.tools.types import check_instance, first_not_None, to_tuple
 from hysop.tools.decorators import debug
-from hysop.fields.continuous_field import Field
+from hysop.fields.continuous_field import ScalarField, Field
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.backend.device.opencl.operator.custom_symbolic import OpenClCustomSymbolicOperator
 from hysop.symbolic.relational import Assignment
@@ -23,7 +23,7 @@ class OpenClAnalyticField(OpenClCustomSymbolicOperator):
 
         Parameters
         ----------
-        field: hysop.field.continuous_field.Field
+        field: hysop.field.continuous_field.ScalarField
             Continuous field to be modified.
         formula : sm.Basic or array-like of sm.Basic
             field.nb_components symbolic expressions as a tuple.
@@ -33,7 +33,7 @@ class OpenClAnalyticField(OpenClCustomSymbolicOperator):
             Base class arguments.
         """
         formula = to_tuple(formula)
-        check_instance(field, Field)
+        check_instance(field, ScalarField)
         check_instance(formula, tuple, values=(type(None),sm.Basic), size=field.nb_components)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
         
diff --git a/hysop/backend/host/python/operator/analytic.py b/hysop/backend/host/python/operator/analytic.py
index e2969cfe0..7d9dd1a93 100644
--- a/hysop/backend/host/python/operator/analytic.py
+++ b/hysop/backend/host/python/operator/analytic.py
@@ -3,7 +3,7 @@ from hysop.tools.types import check_instance, first_not_None
 from hysop.tools.decorators import debug
 from hysop.backend.host.host_operator import HostOperator
 from hysop.core.graph.graph import op_apply
-from hysop.fields.continuous_field import Field
+from hysop.fields.continuous_field import Field, ScalarField
 from hysop.parameters.parameter import Parameter
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 
@@ -24,7 +24,7 @@ class PythonAnalyticField(HostOperator):
 
         Parameters
         ----------
-        field: hysop.field.continuous_field.Field
+        field: hysop.field.continuous_field.ScalarField
             Continuous field to be modified.
         formula : callable
             The formula to be applied onto the field.
@@ -44,20 +44,20 @@ class PythonAnalyticField(HostOperator):
         """
         extra_input_kwds = first_not_None(extra_input_kwds, {})
 
-        check_instance(field, Field)
+        check_instance(field, ScalarField)
         assert callable(formula), type(formula)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
         check_instance(extra_input_kwds, dict, keys=str)
 
         input_fields  = {}
-        output_fields = { field: variables[field] }
+        output_fields = { field: self.get_topo_descriptor(variables, field) }
         input_params  = {}
 
         extra_kwds = {}
         map_fields = {}
         for (k,v) in extra_input_kwds.iteritems():
             if isinstance(v, Field):
-                input_fields[v] = variables[v]
+                input_fields[v] = self.get_topo_descriptor(variables, v)
                 map_fields[v] = k
             elif isinstance(v, Parameter):
                 input_params[k] = v
@@ -68,6 +68,7 @@ class PythonAnalyticField(HostOperator):
         super(PythonAnalyticField, self).__init__(input_fields=input_fields, 
                 output_fields=output_fields,
                 input_params=input_params, **kwds)
+
         self.field = field
         self.formula = formula
         self.extra_kwds = extra_kwds
@@ -83,8 +84,8 @@ class PythonAnalyticField(HostOperator):
         map_fields = self.map_fields
         assert 'data'   not in extra_kwds
         assert 'coords' not in extra_kwds
-        extra_kwds['data']   = dfield.compute_data
-        extra_kwds['coords'] = dfield.get_attributes('compute_mesh_coords')
+        extra_kwds['data']   = dfield.compute_data[0]
+        extra_kwds['coords'] = dfield.compute_mesh_coords
         for (field, dfield) in self.input_discrete_fields.iteritems():
             assert field.name not in extra_kwds, field.name
             extra_kwds[map_fields[field]] = dfield.compute_data
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index c5f969119..9501cf429 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -9,7 +9,7 @@ from hysop.constants import BoundaryCondition
 from hysop.backend.host.host_operator import ComputeGranularity
 from hysop.backend.host.host_directional_operator import HostDirectionalOperator
 from hysop.operator.base.advection_dir import DirectionalAdvectionBase
-from hysop.methods import Interpolation
+from hysop.methods import Interpolation, PolynomialInterpolation
 from hysop.numerics.odesolvers.runge_kutta import ExplicitRungeKutta, Euler, RK2, RK3, RK4
 
 DEBUG = False
@@ -187,7 +187,6 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
         lidx: int32 empty input buffer, destroyed
         ridx: int32 empty input buffer, destroyed
         """
-        assert self.interp is Interpolation.LINEAR # legacy linear interpolator
         N = dX.shape[-1]
 
         Iy = I[:-1]
@@ -431,7 +430,9 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
         super(PythonDirectionalAdvection, self).handle_method(method)
         cr = method.pop(ComputeGranularity)
         assert 0 <= cr <= self.velocity.dim-1
-        assert self.interp == Interpolation.LINEAR
+        msg='Interpolation {}.{} is not supported for operator {}.'.format(
+                self.interp.__class__.__name__, self.interp, self.__class__.__name__)
+        assert self.interp in (Interpolation.LINEAR, PolynomialInterpolation.LINEAR), msg
         self.compute_granularity = cr
 
     @classmethod
diff --git a/hysop/core/graph/computational_node.py b/hysop/core/graph/computational_node.py
index e7b0c86b5..c20dba7c9 100644
--- a/hysop/core/graph/computational_node.py
+++ b/hysop/core/graph/computational_node.py
@@ -18,6 +18,7 @@ from hysop.core.graph.continuous import OperatorBase
 from hysop.topology.topology import Topology, TopologyView
 from hysop.tools.decorators import debug
 from hysop.tools.warning import HysopWarning
+from hysop.topology.cartesian_descriptor import get_topo_descriptor_discretization
      
 
 def base_initialized(f):
@@ -686,6 +687,11 @@ class ComputationalGraphNode(OperatorBase):
         msg=msg.format(field.short_description())
         raise KeyError(msg)
     
+    @classmethod
+    def get_topo_discretization(cls, variables, field):
+        topo = cls.get_topo_descriptor(variables=variables, field=field)
+        return get_topo_descriptor_discretization(topo)
+    
     @classmethod
     def supports_multiple_topologies(cls):
         """
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index 7330ae255..5749e83f8 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -28,7 +28,7 @@ from hysop.core.mpi.topo_tools import TopoTools
 
 class CartesianDiscreteScalarFieldViewContainerI(object):
     def initialize(self, formula, vectorize=False,
-            without_ghosts=False, exchange_ghosts=True,
+            without_ghosts=True, exchange_ghosts=True,
             exchange_kwds=None, only_finite=True,
             reorder=None, quiet=False, components=None, **kwds):
         """
@@ -43,22 +43,23 @@ class CartesianDiscreteScalarFieldViewContainerI(object):
         vectorize : bool, optional
             True if formula must be vectorized
             (i.e. is of type 'user_func_2', see :ref:`fields` for details)
-        only_finite: bool, optional
-            Check that initialized values are not +inf, -inf or NaN.
-            Defaults to True.
+        without_ghosts: boolean, optional, defaults to False
+            Do not initialize ghosts using formula.
+            In this case, ghosts may be initialized by ghost exchange (see next parameter).
         exchange_ghosts: bool, optional, defaults to True
             Should we exchange ghosts after initialization ?
-        without_ghosts: boolean, optional, defaults to False
-            Do not initialize ghosts (only for formula init).
         exchange_kwds: dict, optional,
             Extra exchange keyword arguments passed to ghost exchange.
             Only used if exchange_ghosts is set to True.
+        only_finite: bool, optional
+            Check that initialized values are not +inf, -inf or NaN.
+            Defaults to True.
         reorder: str or tuple of str, optional
             Reorder all kwd according to current topology state.
             ie. kwd should be contained in kwds and kwds[kwd] should
                 be of length self.domain.dim.
         kwds : dict
-            Extra keyword arguments passed to formula, optional, depends on formula.
+            Extra keyword arguments passed to formula, optional, depends on passed formula.
         """
         if (not self.has_unique_backend()):
             msg='Cannot initialize a tensor field defined on multiple backends at a time.'
@@ -185,7 +186,7 @@ class CartesianDiscreteScalarFieldViewContainerI(object):
                 if (ddst is None):
                     continue
                 src_slices = all_src_slices[i]
-                dst_slices = kwds.pop('dst_slices', dfields[i].compute_slices)
+                dst_slices = kwds.pop('dst_slices', Ellipsis)
                 check_instance(dsrc, (np.ndarray, HostArray))
                 check_instance(ddst, (np.ndarray, HostArray))
                 src_shape = dsrc[src_slices].shape
diff --git a/hysop/fields/tests/test_cartesian.py b/hysop/fields/tests/test_cartesian.py
index 6bac6d6af..7b178728b 100644
--- a/hysop/fields/tests/test_cartesian.py
+++ b/hysop/fields/tests/test_cartesian.py
@@ -12,27 +12,23 @@ from hysop.topology.cartesian_topology import CartesianTopology, CartesianTopolo
 from hysop.testsenv import iter_clenv, test_context, domain_boundary_iterator
 from hysop.tools.numerics import is_fp, is_integer
 
-def __random_init(data, coords):
-    shape = data[0].shape
-    dtype = data[0].dtype
+def __random_init(data, coords, component):
+    shape = data.shape
+    dtype = data.dtype
     if is_integer(dtype):
-        for d in data:
-            d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+        data[...] = npw.random.random_integers(low=0, high=255, size=shape)
     elif is_fp(dtype):
-        for d in data:
-            d[...] = npw.random.random(size=d.shape)
+        data[...] = npw.random.random(size=shape)
     else:
         msg = 'Unknown dtype {}.'.format(dtype)
         raise NotImplementedError(msg)
 
-def __zero_init(data, coords):
-    for d in data:
-        d[...] = 0
+def __zero_init(data, coords, component):
+    data[...] = 0
 
 def __cst_init(cst):
-    def __init(data,coords):
-        for d in data:
-            d[...] = cst
+    def __init(data,coords,component):
+        data[...] = cst
     return __init
 
 def test_serial_initialization_1d():
diff --git a/hysop/methods.py b/hysop/methods.py
index c122b82af..741c3f020 100644
--- a/hysop/methods.py
+++ b/hysop/methods.py
@@ -33,7 +33,7 @@ from hysop.constants import Backend, Precision, BoundaryCondition, \
 
 from hysop.operator.spatial_filtering import FilteringMethod
 from hysop.numerics.interpolation.interpolation import Interpolation, MultiScaleInterpolation
-from hysop.numerics.interpolation.polynomial import PolynomialInterpolator
+from hysop.numerics.interpolation.polynomial import PolynomialInterpolator, PolynomialInterpolation
 from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator
 from hysop.numerics.remesh.remesh import Remesh
 from hysop.numerics.splitting.strang import StrangOrder
diff --git a/hysop/numerics/tests/test_fft.py b/hysop/numerics/tests/test_fft.py
index c08d154d5..763f455bf 100644
--- a/hysop/numerics/tests/test_fft.py
+++ b/hysop/numerics/tests/test_fft.py
@@ -19,9 +19,10 @@ from hysop.numerics.fft.fft  import mk_shape, HysopFFTDataLayoutError
 from hysop.numerics.fft.numpy_fft  import NumpyFFT
 from hysop.numerics.fft.scipy_fft  import ScipyFFT
 from hysop.numerics.fft.fftw_fft   import FftwFFT
-from hysop.numerics.fft._mkl_fft   import MklFFT
 from hysop.numerics.fft.gpyfft_fft import GpyFFT
 
+# For now Intel has a bug in their MKL FFTW-like interface so we do not test it
+# from hysop.numerics.fft._mkl_fft import MklFFT
 
 class TestFFT(object):
 
@@ -29,9 +30,9 @@ class TestFFT(object):
         Implementation.PYTHON: {
             'numpy': NumpyFFT(warn_on_allocation=False),
             'scipy': ScipyFFT(warn_on_allocation=False),
-            'mkl':   MklFFT(warn_on_allocation=False),
             'fftw':  FftwFFT(warn_on_allocation=False,
                              warn_on_misalignment=False)
+            #'mkl':   MklFFT(warn_on_allocation=False),
         },
         Implementation.OPENCL: {}
     }
diff --git a/hysop/operator/analytic.py b/hysop/operator/analytic.py
index d1f4a38ce..cff14829f 100644
--- a/hysop/operator/analytic.py
+++ b/hysop/operator/analytic.py
@@ -2,12 +2,105 @@
 """
 from hysop.constants import Backend, Implementation
 from hysop.fields.continuous_field import Field
-from hysop.tools.types import check_instance
+from hysop.tools.types import check_instance, first_not_None, to_tuple
 from hysop.tools.decorators import debug
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
+from hysop.core.graph.node_generator import ComputationalGraphNodeGenerator
 
-class AnalyticField(ComputationalGraphNodeFrontend):
+class AnalyticField(ComputationalGraphNodeGenerator):
+    """
+    Applies an analytic formula, given by user, on all contained scalar fields.
+    Formula may be given in different formats depending on the
+    chosen implementation backend.
+    """
+
+    @debug
+    def __init__(self, field, formula, variables, extra_input_kwds=None,
+            implementation=None, base_kwds=None, **kwds):
+        """
+        AnalyticField operator frontend.
+
+        Apply a user-defined formula onto a field, possibly 
+        dependent on space variables and external fields/parameters.
+
+        Parameters
+        ----------
+        field: hysop.field.continuous_field.Field
+            Continuous field to be modified.
+        formula : python function or sympy expression or tupe
+            The formula to be applied onto the scalar fields.
+            If the formula is a tuple of the length of the number of scalar fields,
+            fomula[component] is used instead.
+        variables: dict
+            Dictionary of fields as keys and topology descriptors as values.
+        implementation: Implementation, optional, defaults to None
+            target implementation, should be contained in available_implementations().
+            If None, implementation will be set to default_implementation().
+        extra_input_kwds: dict, optional
+            Extra inputs that will be forwarded to the formula.
+            Fields and Parameters are handled correctly as input requirements.
+            Only used for Implementation.PYTHON, discarded for other implementations.
+        base_kwds: dict, optional
+            Base class keyword arguments.
+        kwds: dict, optional
+            Extra parameters passed towards operator implementation.
+        """
+        check_instance(field, Field)
+        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+        extra_input_kwds = first_not_None(extra_input_kwds, {})
+        base_kwds = first_not_None(base_kwds, {})
+            
+        assert 'extra_input_kwds' not in kwds
+        assert 'component' not in kwds
+        assert 'coords' not in kwds
+
+        if (implementation is Implementation.PYTHON) and (extra_input_kwds is not None):
+            candidate_input_tensors = filter(lambda f: isinstance(f, Field), extra_input_kwds.values())
+        else:
+            extra_input_kwds = {}
+            candidate_input_tensors = ()
+        candidate_output_tensors = (field,)
+        
+        formula = to_tuple(formula)
+        if len(formula) == 1:
+            formula = formula*field.nb_components
+        check_instance(formula, tuple, size=field.nb_components)
+        
+        super(AnalyticField, self).__init__(
+                candidate_input_tensors=candidate_input_tensors,
+                candidate_output_tensors=candidate_output_tensors,
+                **base_kwds)
+
+        self._fields    = field.fields
+        self._formula   = formula
+        self._variables = variables
+        self._extra_input_kwds = extra_input_kwds
+        self._implementation = implementation
+        self._kwds = kwds
+    
+    @debug
+    def _generate(self):
+        nodes = []
+        impl      = self._implementation
+        variables = self._variables
+        assert len(self._formula)==len(self._fields)
+        for component, (formula, field) in enumerate(zip(self._formula, self._fields)):
+            kwds = self._kwds.copy()
+            extra_input_kwds = self._extra_input_kwds.copy()
+            extra_input_kwds['component'] = component
+            node = AnalyticScalarField(
+                    field=field,
+                    formula=formula,
+                    variables=variables,
+                    implementation=impl, 
+                    extra_input_kwds=extra_input_kwds,
+                    **kwds)
+            nodes.append(node)
+        return nodes
+
+
+class AnalyticScalarField(ComputationalGraphNodeFrontend):
     """
     Applies an analytic formula, given by user, on its field.
     Formula may be given in different formats depending on the
@@ -63,6 +156,6 @@ class AnalyticField(ComputationalGraphNodeFrontend):
         if (implementation is Implementation.PYTHON):
             kwds['extra_input_kwds'] = extra_input_kwds
 
-        super(AnalyticField, self).__init__(field=field, formula=formula, 
+        super(AnalyticScalarField, self).__init__(field=field, formula=formula, 
                 variables=variables, implementation=implementation, 
                 base_kwds=base_kwds, **kwds)
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index 30db42947..cb23b75f5 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -332,7 +332,6 @@ class RemeshRestrictionFilterBase(RestrictionFilterBase):
         assert abs(weights.sum() - 1.0) < 1e-8, weights.sum()
         assert abs(npw.sum(nz_weights.values()) - 1.0) < 1e-8, npw.sum(nz_weights.values()) 
 
-        self.iratio     = iratio
         self.weights    = weights
         self.nz_weights = nz_weights
 
@@ -342,12 +341,12 @@ class RemeshRestrictionFilterBase(RestrictionFilterBase):
             return
         super(RemeshRestrictionFilterBase, self).discretize()
         dFin, dFout  = self.dFin, self.dFout
-
-        iratio =  dFin.compute_resolution / dFout.compute_resolution
-        self.compute_weights(iratio)
+    
+        grid_ratio = self.grid_ratio
+        self.compute_weights(grid_ratio)
         
         remesh_ghosts    = self.remesh_ghosts
-        fine_grid_ghosts = iratio*remesh_ghosts - 1
+        fine_grid_ghosts = grid_ratio*remesh_ghosts - 1
         fin  = dFin.sdata[dFin.local_slices(ghosts=fine_grid_ghosts)]
         fout = dFout.compute_buffers[0]
 
@@ -368,7 +367,7 @@ class SubgridRestrictionFilterBase(RestrictionFilterBase):
         dFin, dFout  = self.dFin, self.dFout
 
         grid_ratio = self.grid_ratio
-        view = tuple(slice(None,None,r) for r in iratio)
+        view = tuple(slice(None,None,r) for r in grid_ratio)
         
         fin  = dFin.compute_buffers[0][view]
         fout = dFout.compute_buffers[0]
@@ -376,7 +375,7 @@ class SubgridRestrictionFilterBase(RestrictionFilterBase):
         msg='Something went wrong during slicing: fin.shape={}, fout.shape={}'
         msg=msg.format(fin.shape, fout.shape)
         assert (fin.shape == fout.shape), msg
-        assert npw.prod(iratio) == npw.prod(self.iratio), msg
+        assert npw.prod(grid_ratio) == npw.prod(self.iratio), msg
 
         self.fin, self.fout = fin, fout
 
diff --git a/hysop/operator/spatial_filtering.py b/hysop/operator/spatial_filtering.py
index f485aec69..7facb047c 100644
--- a/hysop/operator/spatial_filtering.py
+++ b/hysop/operator/spatial_filtering.py
@@ -22,7 +22,7 @@ class SpatialFilterFrontend(MultiComputationalGraphNodeFrontend):
                  base_kwds=None,
                  **kwds):
         """
-        Initialize a RestrictionFilter operator.
+        Initialize a SpatialFilter operator.
 
         Parameters
         ----------
@@ -43,7 +43,7 @@ class SpatialFilterFrontend(MultiComputationalGraphNodeFrontend):
 
         Notes
         -----
-        An implementation should at least support the RestrictionFilterBase interface.
+        An implementation should at least support the SpatialFilterBase interface.
         """
         check_instance(input_variable, tuple, size=2)
         check_instance(output_variable, tuple, size=2)
@@ -209,16 +209,15 @@ class SpatialFilter(ComputationalGraphNodeGenerator):
     def _generate(self):
         nodes = []
         for (ifield, ofield) in zip(self._input_fields, self._output_fields):
-            stopo = ComputationalGraphNode.get_topo_descriptor(self._input_variables, ifield)
-            ttopo = ComputationalGraphNode.get_topo_descriptor(self._output_variables, ofield)
-            fm    = self._fm
-            impl  = self._impl
-            kwds  = self._kwds.copy()
-            
-            # TODO handle other topology descriptors
+            stopo = ComputationalGraphNode.get_topo_discretization(self._input_variables, ifield)
+            ttopo = ComputationalGraphNode.get_topo_discretization(self._output_variables, ofield)
             check_instance(stopo, tuple, values=(int,long))
             check_instance(ttopo, tuple, values=(int,long))
             assert len(stopo)==len(ttopo)
+
+            fm    = self._fm
+            impl  = self._impl
+            kwds  = self._kwds.copy()
             
             # if source topology is destination topology there is nothing to be done
             if (ttopo == stopo):
diff --git a/hysop/operator/tests/test_absorption.py b/hysop/operator/tests/test_absorption.py
index 8da2141d2..c3a8354de 100644
--- a/hysop/operator/tests/test_absorption.py
+++ b/hysop/operator/tests/test_absorption.py
@@ -47,20 +47,17 @@ class TestVorticityAbsorption(object):
         pass
 
     @staticmethod
-    def __random_init(data, coords):
-        dtype = data[0].dtype
+    def __random_init(data, coords, component):
+        dtype = data.dtype
         if is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape).astype(dtype=dtype)
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
     @staticmethod
-    def __velo_init(data, coords):
-        data[0][...] = 1.
-        data[1][...] = 0.
-        data[2][...] = 0.
+    def __velo_init(data, coords, component):
+        data[...] = [1,0,0][component]
 
     def _test(self, dim, dtype,
               size_min=None, size_max=None):
diff --git a/hysop/operator/tests/test_analytic.py b/hysop/operator/tests/test_analytic.py
index 47b66cd7e..83e397047 100644
--- a/hysop/operator/tests/test_analytic.py
+++ b/hysop/operator/tests/test_analytic.py
@@ -12,7 +12,7 @@ from hysop.tools.types import check_instance, first_not_None
 from hysop.tools.io_utils import IO
 from hysop.tools.numpywrappers import npw
 from hysop.parameters.scalar_parameter import ScalarParameter
-from hysop.operator.analytic import AnalyticField, Implementation
+from hysop.operator.analytic import AnalyticField, AnalyticScalarField, Implementation
 
 from hysop import Field, Box
 
@@ -87,20 +87,18 @@ class TestAnalyticField(object):
         pass
 
     @staticmethod
-    def __random_init(data, coords):
-        dtype = data[0].dtype
+    def __random_init(data, coords, component):
+        dtype = data.dtype
         if is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape).astype(dtype=dtype)
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
-    @classmethod
-    def __analytic_init(cls, data, coords, fns, t):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            d[...] = npw.asarray(fn(*(coord+(t(),)))).astype(d.dtype)
+    @staticmethod
+    def __analytic_init(data, coords, fns, t, component):
+        fn = fns[component]
+        data[...] = npw.asarray(fn(*(coords+(t(),)))).astype(data.dtype)
 
     def _test(self, dim, dtype,
             size_min=None, size_max=None):
@@ -131,7 +129,7 @@ class TestAnalyticField(object):
         print ' >Parameter t has been set to {}.'.format(self.t())
         print ' >Testing all implementations:'
 
-        implementations = AnalyticField.implementations()
+        implementations = AnalyticScalarField.implementations()
 
         variables = { F:shape }
         fns   = self.analytic_functions[dim]['F']
@@ -177,7 +175,7 @@ class TestAnalyticField(object):
                 dF.initialize(self.__random_init)
                 op.apply()
 
-                Fout   = tuple( data.get().handle.copy() for data in dF.data )
+                Fout = tuple( data.get().handle.copy() for data in dF.data )
                 self._check_output(impl, op, Fref, Fout)
 
     @classmethod
@@ -215,7 +213,7 @@ class TestAnalyticField(object):
 
             print
             print
-            print 'Test output comparisson for {} failed for component {}:'.format(name, i)
+            print 'Test output comparisson for {} failed for component {}:'.format(iname, i)
             print ' *has_nan: {}'.format(has_nan)
             print ' *has_inf: {}'.format(has_inf)
             print ' *dinf={} ({} eps)'.format(dinf, deps)
@@ -247,7 +245,7 @@ class TestAnalyticField(object):
                         print
                 print
 
-            msg = 'Test failed for {} on component {} for implementation {}.'.format(name, i, impl)
+            msg = 'Test failed for {} on component {} for implementation {}.'.format(iname, i, impl)
             raise RuntimeError(msg)
 
 
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
index a0e95366a..d326ddbb3 100644
--- a/hysop/operator/tests/test_diffusion.py
+++ b/hysop/operator/tests/test_diffusion.py
@@ -60,26 +60,24 @@ class TestDiffusionOperator(object):
             self._test_one(shape=shape, dim=dim, dtype=dtype,
                     is_inplace=is_inplace, domain=domain,
                     Fin=Fin, Fout=Fout, nu=nu)
-
+    
     @staticmethod
-    def __random_init(data, coords, dtype):
+    def __random_init(data, coords, dtype, component):
+        shape = data.shape
         if is_integer(dtype):
-            for d in data:
-                d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+            data[...] = npw.random.random_integers(low=0, high=255, size=shape)
         elif is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape)
+            data[...] = npw.random.random(size=data.shape)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
     @staticmethod
-    def __scalar_init(data, coords, dtype):
+    def __scalar_init(data, coords, dtype, component):
         if is_fp(dtype):
-            for i,(d,c) in enumerate(zip(data,coords)):
-                d[...] = 1
-                for xi in c:
-                    d[...] *= npw.cos(xi*(i+1))
+            data[...] = 1
+            for xi in coords:
+                data[...] *= npw.cos(xi*(component+1))
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
diff --git a/hysop/operator/tests/test_directional_advection.py b/hysop/operator/tests/test_directional_advection.py
index b424ab569..c1e4f89bb 100644
--- a/hysop/operator/tests/test_directional_advection.py
+++ b/hysop/operator/tests/test_directional_advection.py
@@ -236,7 +236,7 @@ class TestDirectionalAdvectionOperator(object):
                                 for i in xrange(dsout.nb_components):
                                     di = npw.abs(reference[i] - output[i])
                                     max_di = npw.max(di)
-                                    neps = 500
+                                    neps = 1000
                                     max_tol = neps*npw.finfo(dsout.dtype).eps
                                     if (max_di>max_tol):
                                         print 'FATAL ERROR: Could not match other implementation results.'
diff --git a/hysop/operator/tests/test_directional_diffusion.py b/hysop/operator/tests/test_directional_diffusion.py
index 4949daad9..626192604 100644
--- a/hysop/operator/tests/test_directional_diffusion.py
+++ b/hysop/operator/tests/test_directional_diffusion.py
@@ -77,25 +77,22 @@ class TestDirectionalDiffusionOperator(object):
                             Fin=Fin, Fout=Fout, coeffs=coeffs)
 
     @staticmethod
-    def __random_init(data, coords, dtype):
-        shape = data[0].shape
+    def __random_init(data, coords, dtype, component):
+        shape = data.shape
         if is_integer(dtype):
-            for d in data:
-                d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+            data[...] = npw.random.random_integers(low=0, high=255, size=shape)
         elif is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape)
+            data[...] = npw.random.random(size=data.shape)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
     @staticmethod
-    def __scalar_init(data, coords, dtype):
+    def __scalar_init(data, coords, dtype, component):
         if is_fp(dtype):
-            for i,(d,coord) in enumerate(zip(data,coords)):
-                d[...] = 1
-                for xi in coord:
-                    d[...] *= npw.cos(xi*(i+1))
+            data[...] = 1
+            for xi in coords:
+                data[...] *= npw.cos(xi*(component+1))
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
diff --git a/hysop/operator/tests/test_directional_stretching.py b/hysop/operator/tests/test_directional_stretching.py
index 326327d3a..48e40dcf9 100644
--- a/hysop/operator/tests/test_directional_stretching.py
+++ b/hysop/operator/tests/test_directional_stretching.py
@@ -85,15 +85,13 @@ class TestDirectionalStretchingOperator(object):
                                 C=C, A=A, formulation=formulation)
 
     @staticmethod
-    def __random_init(data, coords):
-        shape = data[0].shape
-        dtype = data[0].dtype
+    def __random_init(data, coords, component):
+        shape = data.shape
+        dtype = data.dtype
         if is_integer(dtype):
-            for d in data:
-                d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+            data[...] = npw.random.random_integers(low=0, high=255, size=shape)
         elif is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape)
+            data[...] = npw.random.random(size=shape)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
diff --git a/hysop/operator/tests/test_fd_derivative.py b/hysop/operator/tests/test_fd_derivative.py
index 6bca7fbe3..12f7f188d 100644
--- a/hysop/operator/tests/test_fd_derivative.py
+++ b/hysop/operator/tests/test_fd_derivative.py
@@ -93,20 +93,17 @@ class TestFiniteDifferencesDerivative(object):
         pass
 
     @staticmethod
-    def __random_init(data, coords):
-        dtype = data[0].dtype
+    def __random_init(data, coords, component):
+        dtype = data.dtype
         if is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape).astype(dtype=dtype)
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
     @classmethod
-    def __analytic_init(cls, data, coords, fns, t):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            d[...] = fn(*(coord+(t(),))).astype(d.dtype)
+    def __analytic_init(cls, data, coords, component, fns, t):
+        data[...] = fns[component](*(coords+(t(),))).astype(data.dtype)
 
     def _test(self, dim, dtype,
             size_min=None, size_max=None):
diff --git a/hysop/operator/tests/test_poisson.py b/hysop/operator/tests/test_poisson.py
index 4a5a04d71..693ed2040 100644
--- a/hysop/operator/tests/test_poisson.py
+++ b/hysop/operator/tests/test_poisson.py
@@ -75,21 +75,17 @@ class TestPoissonOperator(object):
         return (analytic_expressions, analytic_functions)
 
     @staticmethod
-    def __random_init(data, coords, dtype):
-        for d in data:
-            if is_fp(d.dtype):
-                d[...] = npw.random.random(size=d.shape).astype(dtype=d.dtype)
-            else:
-                msg = 'Unknown dtype {}.'.format(d.dtype)
-                raise NotImplementedError(msg)
+    def __random_init(data, coords, component, dtype):
+        if is_fp(dtype):
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
+        else:
+            msg = 'Unknown dtype {}.'.format(dtype)
+            raise NotImplementedError(msg)
 
     @staticmethod
-    def __analytic_init(data, coords, dtype, fns):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            coord = tuple(c.astype(d.dtype) for c in coord)
-            d[...] = fn(*coord).astype(d.dtype)
-
+    def __analytic_init(data, coords, fns, component, dtype):
+        fn = fns[component]
+        data[...] = fn(*coords).astype(data.dtype)
 
     def _test(self, dim, dtype, max_runs=5,
                    polynomial=False, size_min=None, size_max=None):
diff --git a/hysop/operator/tests/test_poisson_curl.py b/hysop/operator/tests/test_poisson_curl.py
index 82d7bdd3d..30ec59dae 100644
--- a/hysop/operator/tests/test_poisson_curl.py
+++ b/hysop/operator/tests/test_poisson_curl.py
@@ -84,22 +84,17 @@ class TestPoissonCurlOperator(object):
 
 
     @staticmethod
-    def __random_init(data, coords, dtype):
-        for d in data:
-            if is_fp(d.dtype):
-                d[...] = npw.random.random(size=d.shape).astype(dtype=d.dtype)
-            else:
-                msg = 'Unknown dtype {}.'.format(d.dtype)
-                raise NotImplementedError(msg)
+    def __random_init(data, coords, component, dtype):
+        if is_fp(dtype):
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
+        else:
+            msg = 'Unknown dtype {}.'.format(dtype)
+            raise NotImplementedError(msg)
 
     @staticmethod
-    def __analytic_init(data, coords, dtype, fns):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            coord = tuple(c.astype(d.dtype) for c in coord)
-            d[...] = fn(*coord).astype(d.dtype)
-
-
+    def __analytic_init(data, coords, fns, component, dtype):
+        fn = fns[component]
+        data[...] = fn(*coords).astype(data.dtype)
 
     def _test(self, dim, dtype, max_runs=5,
             polynomial=False, size_min=None, size_max=None):
diff --git a/hysop/operator/tests/test_restriction_filter.py b/hysop/operator/tests/test_restriction_filter.py
index dcf464c3b..744351cc9 100755
--- a/hysop/operator/tests/test_restriction_filter.py
+++ b/hysop/operator/tests/test_restriction_filter.py
@@ -5,7 +5,7 @@ from hysop.testsenv import __ENABLE_LONG_TESTS__
 from hysop.tools.io_utils import IO
 from hysop.tools.numpywrappers import npw
 from hysop.tools.types import first_not_None
-from hysop.operator.spatial_filtering import RestrictionFilter
+from hysop.operator.spatial_filtering import SpatialFilter
 from hysop.methods import FilteringMethod
 from hysop.topology.cartesian_topology import CartesianTopology
 from hysop.constants import implementation_to_backend, Implementation, HYSOP_REAL
@@ -18,17 +18,16 @@ from hysop import Field, Box, MPIParams
 class TestMultiresolutionFilter(object):
 
     @staticmethod
-    def __f_init(data, coords):
+    def __f_init(data, coords, component):
         from numpy import sin, cos
-        (x, y, z) = coords[0]
-        data[0][...] = - cos(x) * sin(y) * sin(z)
+        (x, y, z) = coords
+        data[...] = - cos(x) * sin(y) * sin(z)
 
     @staticmethod
-    def __random_init(data, coords):
-        dtype = data[0].dtype
+    def __random_init(data, coords, component):
+        dtype = data.dtype
         if is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape).astype(dtype=dtype)
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
@@ -93,7 +92,9 @@ class TestMultiresolutionFilter(object):
             if impl is Implementation.PYTHON:
                 msg='   *Python: '
                 print msg,
-                yield RestrictionFilter(input_variables={f: topo_f}, output_variables={f: topo_c}, **base_kwds)
+                yield SpatialFilter(input_variables={f: topo_f},
+                        output_variables={f: topo_c}, 
+                        **base_kwds)
                 print
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
diff --git a/hysop/operator/tests/test_solenoidal_projection.py b/hysop/operator/tests/test_solenoidal_projection.py
index 234a1e845..e9eed1747 100644
--- a/hysop/operator/tests/test_solenoidal_projection.py
+++ b/hysop/operator/tests/test_solenoidal_projection.py
@@ -149,28 +149,24 @@ class TestSolenoidalProjectionOperator(object):
             print '>> TESTED ALL {}D BOUNDARY CONDITIONS <<'.format(dim)
             print
             print
-
-    @classmethod
-    def __analytic_init(cls, data, coords, dtype, fns):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            d[...] = npw.asarray(fn(*coord)).astype(dtype)
-
-    @classmethod
-    def __zero_init(cls, data, coords, dtype):
-        for d in data:
-            d[...] = 0
-
+    
     @staticmethod
-    def __random_init(data, coords, dtype):
-        shape = data[0].shape
+    def __random_init(data, coords, component, dtype):
         if is_fp(dtype):
-            for d in data:
-                d[...] = npw.random.random(size=d.shape).astype(dtype=dtype)
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
         else:
             msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
+    @staticmethod
+    def __analytic_init(data, coords, fns, component, dtype):
+        fn = fns[component]
+        data[...] = npw.asarray(fn(*coords)).astype(dtype)
+
+    @classmethod
+    def __zero_init(cls, data, coords, dtype, component):
+        data[...] = 0
+
     def _test_one(self, shape, dtype, polynomial,
             domain, U, U0, U1, divU, divU0, divU1):
 
diff --git a/hysop/operator/tests/test_spectral_curl.py b/hysop/operator/tests/test_spectral_curl.py
index 3407397bd..273719db8 100644
--- a/hysop/operator/tests/test_spectral_curl.py
+++ b/hysop/operator/tests/test_spectral_curl.py
@@ -79,24 +79,18 @@ class TestSpectralCurl(object):
         analytic_functions   = {'Fin':fFins, 'Fout':fFouts}
         return (analytic_expressions, analytic_functions)
 
-    
-    @staticmethod
-    def __random_init(data, coords, dtype):
-        for d in data:
-            if is_fp(d.dtype):
-                d[...] = npw.random.random(size=d.shape).astype(dtype=d.dtype)
-            else:
-                msg = 'Unknown dtype {}.'.format(d.dtype)
-                raise NotImplementedError(msg)
-    
     @staticmethod
-    def __analytic_init(data, coords, dtype, fns):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            coord = tuple(c.astype(d.dtype) for c in coord)
-            d[...] = fn(*coord).astype(d.dtype)
-
+    def __random_init(data, coords, component, dtype):
+        if is_fp(dtype):
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
+        else:
+            msg = 'Unknown dtype {}.'.format(dtype)
+            raise NotImplementedError(msg)
 
+    @staticmethod
+    def __analytic_init(data, coords, fns, component, dtype):
+        fn = fns[component]
+        data[...] = npw.asarray(fn(*coords)).astype(dtype)
 
     def _test(self, dim, dtype, nb_components, max_runs=5,
             polynomial=False, size_min=None, size_max=None):
diff --git a/hysop/operator/tests/test_spectral_derivative.py b/hysop/operator/tests/test_spectral_derivative.py
index 6a5b2b96e..b8f3cb7f0 100644
--- a/hysop/operator/tests/test_spectral_derivative.py
+++ b/hysop/operator/tests/test_spectral_derivative.py
@@ -92,20 +92,18 @@ class TestSpectralDerivative(object):
         pass
 
     @staticmethod
-    def __random_init(data, coords):
-        for d in data:
-            if is_fp(d.dtype):
-                    d[...] = npw.random.random(size=d.shape).astype(dtype=d.dtype)
-            else:
-                msg = 'Unknown dtype {}.'.format(d.dtype)
-                raise NotImplementedError(msg)
+    def __random_init(data, coords, component):
+        dtype = data.dtype
+        if is_fp(dtype):
+            data[...] = npw.random.random(size=data.shape).astype(dtype=dtype)
+        else:
+            msg = 'Unknown dtype {}.'.format(dtype)
+            raise NotImplementedError(msg)
 
     @staticmethod
-    def __analytic_init(data, coords, fns, t):
-        assert len(fns) == len(data)
-        for (d,fn,coord) in zip(data,fns,coords):
-            coord = tuple(c.astype(d.dtype) for c in coord)
-            d[...] = fn(*(coord+(t(),))).astype(d.dtype)
+    def __analytic_init(data, coords, fns, t, component):
+        fn = fns[component]
+        data[...] = npw.asarray(fn(*(coords+(t(),)))).astype(data.dtype)
 
     def _test(self, dim, dtype, polynomial, max_derivative=2,
             size_min=None, size_max=None, max_runs=None):
diff --git a/hysop/operator/tests/test_transpose.py b/hysop/operator/tests/test_transpose.py
index 6ebe79df3..655f2a73b 100644
--- a/hysop/operator/tests/test_transpose.py
+++ b/hysop/operator/tests/test_transpose.py
@@ -80,19 +80,16 @@ class TestTransposeOperator(object):
                             domain=domain, Fin=Fin, Fout=Fout)
 
     @classmethod
-    def __field_init(cls, data, coords, dtype):
-        shape = data[0].shape
+    def __field_init(cls, data, coords, component, dtype):
+        shape = data.shape
         if is_integer(dtype):
-            for d in data:
-                d[...] = np.random.random_integers(low=0, high=255, size=shape) 
+            data[...] = np.random.random_integers(low=0, high=255, size=shape) 
         elif is_fp(dtype):
-            for d in data:
-                d[...] = np.random.random(size=shape) 
+            data[...] = np.random.random(size=shape) 
         elif is_complex(dtype):
-            for d in data:
-                real = np.random.random(size=shape) 
-                imag = np.random.random(size=shape) 
-                d[...] = real + 1j*imag
+            real = np.random.random(size=shape) 
+            imag = np.random.random(size=shape) 
+            data[...] = real + 1j*imag
         else:
             msg='Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
@@ -223,14 +220,14 @@ class TestTransposeOperator(object):
 
     def perform_tests(self):
         self.test_2d_out_of_place()
+        self.test_3d_out_of_place()
         if __ENABLE_LONG_TESTS__:
-            self.test_3d_out_of_place()
             self.test_4d_out_of_place()
             self.test_upper_dimensions_out_of_place()
         
         self.test_2d_inplace()
+        self.test_3d_inplace()
         if __ENABLE_LONG_TESTS__:
-            self.test_3d_inplace()
             self.test_4d_inplace()
             self.test_upper_dimensions_inplace()
     
diff --git a/hysop/testsenv.py b/hysop/testsenv.py
index 8e1b7e1ec..046f2b4d3 100644
--- a/hysop/testsenv.py
+++ b/hysop/testsenv.py
@@ -93,17 +93,14 @@ if __HAS_OPENCL_BACKEND__:
             cl_environments[cl_device_type] = []
             if (cl_device_type is None):
                 mpi_params = default_mpi_params()
-                if all_platforms:
-                    for i,plat in enumerate(cl.get_platforms()):
-                        for j,dev in enumerate(plat.get_devices()):
-                            cl_env = get_or_create_opencl_env(platform_id=i, device_id=j, 
-                                    mpi_params=mpi_params, **kwds)
+                for i,plat in enumerate(cl.get_platforms()):
+                    for j,dev in enumerate(plat.get_devices()):
+                        cl_env = get_or_create_opencl_env(platform_id=i, device_id=j, 
+                                mpi_params=mpi_params, **kwds)
+                        if (i==__DEFAULT_PLATFORM_ID__) and (j==__DEFAULT_DEVICE_ID__):
+                            cl_environments[None].insert(0, cl_env)
+                        else:
                             cl_environments[None].append(cl_env)
-                else:
-                    cl_env = get_or_create_opencl_env(platform_id=__DEFAULT_PLATFORM_ID__, 
-                            device_id=__DEFAULT_DEVICE_ID__, 
-                            mpi_params=mpi_params, **kwds)
-                    cl_environments[None].append(cl_env)
             else:
                 for cl_env in iter_clenv(cl_device_type=None, all_platforms=True):
                     if (cl_env.device.type & cl_device_type):
@@ -120,6 +117,7 @@ if __HAS_OPENCL_BACKEND__:
             yield cl_env
             if not all_platforms:
                 return
+            
 else:
     opencl_failed = pytest.mark.xfail
     iter_clenv = None
diff --git a/hysop/topology/cartesian_descriptor.py b/hysop/topology/cartesian_descriptor.py
index 80a3315d3..344a44dba 100644
--- a/hysop/topology/cartesian_descriptor.py
+++ b/hysop/topology/cartesian_descriptor.py
@@ -1,5 +1,5 @@
 
-from hysop.tools.types import check_instance
+from hysop.tools.types import check_instance, to_tuple
 from hysop.topology.topology_descriptor import TopologyDescriptor
 from hysop.topology.cartesian_topology import CartesianTopology
 from hysop.tools.parameters import CartesianDiscretization
@@ -209,3 +209,19 @@ Instance of those types can be used to create a CartesianTopologyDescriptor.
 Thus they can be passed in the variables of each operator supporting
 CartesianTopology topologies.
 """
+
+def get_topo_descriptor_discretization(td):
+    """
+    Get grid resolution from any type of CartesianTopologyDescriptor.
+    """
+    check_instance(td, CartesianTopologyDescriptors, allow_none=True)
+    if (td is None):
+        return None
+    elif isinstance(td, CartesianTopology):
+        td = td.grid_resolution
+    elif isinstance(td, CartesianTopologyDescriptor):
+        td = td.grid_resolution
+    elif isinstance(td, CartesianDiscretization):
+        td = td.grid_resolution
+    return to_tuple(td)
+
-- 
GitLab