diff --git a/examples/example_utils.py b/examples/example_utils.py
index f87668a13aca026097a1a35d609efc370d0af744..46cbc096c1ad778eb0201d38e10e52fc41bab144 100644
--- a/examples/example_utils.py
+++ b/examples/example_utils.py
@@ -546,9 +546,13 @@ class HysopArgParser(argparse.ArgumentParser):
                 dest='interpolation_filter',
                 help='Set the default interpolation formula to compute subgrid field values.')
         method.add_argument('-restrict', '--restriction-filter', type=str, 
-                default='subgrid', 
+                default='polynomial', 
                 dest='restriction_filter',
-                help='Set the default filter to pass from fine grids to coarse ones.')
+                help='Set the default restriction formula to pass from fine grids to coarse ones.')
+        method.add_argument('-pi', '--polynomial-interpolator', type=str, 
+                default='LINEAR', 
+                dest='polynomial_interpolator',
+                help='Set the default polynomial interpolator for polynomial interpolation or restriction methods.')
         method.add_argument('-sf', '--stretching-formulation', type=str, 
                 default='conservative', 
                 dest='stretching_formulation',
@@ -567,7 +571,8 @@ class HysopArgParser(argparse.ArgumentParser):
         self._check_default(args, ('compute_granularity', 'fd_order', 'strang_order',
                                     'reprojection_frequency'), int, allow_none=False)
         self._check_default(args, ('time_integrator', 'remesh_kernel', 'interpolation_filter', 
-                                'restriction_filter', 'stretching_formulation'), str, allow_none=False)
+                                'restriction_filter', 'stretching_formulation', 'polynomial_interpolator'), 
+                                str, allow_none=False)
         self._check_positive(args, 'compute_granularity', strict=False, allow_none=False)
         self._check_positive(args, 'fd_order', strict=True, allow_none=False)
         self._check_positive(args, 'strang_order', strict=True, allow_none=False)
@@ -581,10 +586,14 @@ class HysopArgParser(argparse.ArgumentParser):
                 args.remesh_kernel)
         args.stretching_formulation = self._convert_stretching_formulation(
                 'stretching_formulation', args.stretching_formulation)
-        args.interpolation_filter = self._convert_filtering_method('interpolation_filter', 
-                args.interpolation_filter)
-        args.restriction_filter = self._convert_filtering_method('restriction_filter', 
-                args.restriction_filter)
+        args.interpolation_filter = \
+                self._convert_filtering_method('interpolation_filter',  args.interpolation_filter,
+                        allow_subgrid=True)
+        args.restriction_filter = \
+                self._convert_filtering_method('restriction_filter', args.restriction_filter,
+                        allow_subgrid=False)
+        args.polynomial_interpolator = \
+                self._convert_polynomial_interpolation('polynomial_interpolator', args.polynomial_interpolator)
         
     def _add_threading_args(self):
         threading = self.add_argument_group('threading parameters')
@@ -1323,6 +1332,29 @@ class HysopArgParser(argparse.ArgumentParser):
         }
         return self._check_convert(argname, remesh_kernel, remesh_kernels)
     
+    def _convert_polynomial_interpolation(self, argname, pi):
+        from hysop.numerics.interpolation.polynomial import PolynomialInterpolation
+        polynomial_interpolations = {
+          'linear':         PolynomialInterpolation.LINEAR,
+          'cubic':          PolynomialInterpolation.CUBIC,
+          'quintic':        PolynomialInterpolation.QUINTIC,
+          'septic':         PolynomialInterpolation.SEPTIC,
+          'nonic':          PolynomialInterpolation.NONIC,
+          'cubic_fdc2':     PolynomialInterpolation.CUBIC_FDC2,
+          'cubic_fdc4':     PolynomialInterpolation.CUBIC_FDC4,
+          'cubic_fdc6':     PolynomialInterpolation.CUBIC_FDC6,
+          'quintic_fdc2':   PolynomialInterpolation.QUINTIC_FDC2,
+          'quintic_fdc4':   PolynomialInterpolation.QUINTIC_FDC4,
+          'quintic_fdc6':   PolynomialInterpolation.QUINTIC_FDC6,
+          'septic_fdc2':    PolynomialInterpolation.SEPTIC_FDC2,
+          'septic_fdc4':    PolynomialInterpolation.SEPTIC_FDC4,
+          'septic_fdc6':    PolynomialInterpolation.SEPTIC_FDC6,
+          'nonic_fdc2':     PolynomialInterpolation.NONIC_FDC2,
+          'nonic_fdc4':     PolynomialInterpolation.NONIC_FDC4,
+          'nonic_fdc6':     PolynomialInterpolation.NONIC_FDC6,
+        }
+        return self._check_convert(argname, pi, polynomial_interpolations)
+    
     def _convert_stretching_formulation(self, argname, stretching_formulation):
         from hysop.constants import StretchingFormulation
         stretching_formulations = {
@@ -1333,15 +1365,16 @@ class HysopArgParser(argparse.ArgumentParser):
         }
         return self._check_convert(argname, stretching_formulation, stretching_formulations)
     
-    def _convert_filtering_method(self, argname, restriction):
+    def _convert_filtering_method(self, argname, fm, allow_subgrid=True):
         from hysop.methods import FilteringMethod
-        restrictions = {
-            'subgrid':    FilteringMethod.SUBGRID,
+        filtering_methods = {
             'spectral':   FilteringMethod.SPECTRAL,
             'remesh':     FilteringMethod.REMESH,
             'polynomial': FilteringMethod.POLYNOMIAL
         }
-        return self._check_convert(argname, restriction, restrictions)
+        if allow_subgrid:
+            filtering_methods['subgrid'] = FilteringMethod.SUBGRID
+        return self._check_convert(argname, fm, filtering_methods)
         
 class HysopHelpFormatter(ColorHelpFormatter):
     def _format_args(self, *args, **kwds):
diff --git a/examples/multiresolution/scalar_advection.py b/examples/multiresolution/scalar_advection.py
index 9acfdebc216f35f383f8996061f39952e4fd944b..885d756dda11b3a8443e047fa9a31f7526db8b65 100644
--- a/examples/multiresolution/scalar_advection.py
+++ b/examples/multiresolution/scalar_advection.py
@@ -9,7 +9,8 @@ def compute(args):
     from hysop.operators import Advection, DirectionalAdvection, StrangSplitting, \
                                 SpatialFilter, HDF_Writer
     from hysop.methods import Remesh, TimeIntegrator, \
-                              Interpolation, FilteringMethod
+                              Interpolation, FilteringMethod, \
+                              PolynomialInterpolator
 
     ## IO paths
     spectral_path = IO.default_path() + '/spectral'
@@ -46,8 +47,9 @@ def compute(args):
     # Setup operator method dictionnary
     # Advection-Remesh operator discretization parameters
     method = { 
-               TimeIntegrator:      args.time_integrator,
-               Remesh:              args.remesh_kernel,
+       TimeIntegrator:         args.time_integrator,
+       Remesh:                 args.remesh_kernel,
+       PolynomialInterpolator: args.polynomial_interpolator,
     }
     
     # Setup implementation specific variables
diff --git a/hysop/backend/device/opencl/operator/spatial_filtering.py b/hysop/backend/device/opencl/operator/spatial_filtering.py
index d6f168e6e38e4f1646a18cbd23ed39ba6d064ce9..0b4271c79f284de05c133e20172964a2e8046e07 100644
--- a/hysop/backend/device/opencl/operator/spatial_filtering.py
+++ b/hysop/backend/device/opencl/operator/spatial_filtering.py
@@ -51,8 +51,8 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
         for idx in np.ndindex(*gr):
             e = Assignment(fout(gr*I+idx), Fout_values[idx])
             exprs += (e,)
-        
-        interpolate_grid_kernel, _ = ekg.elementwise_kernel('interpolate_grid', 
+        kname='interpolate_grid_{}'.format(self.polynomial_interpolation_method).lower()
+        interpolate_grid_kernel, _ = ekg.elementwise_kernel(kname,
                 *exprs, compute_resolution=self.iter_shape, debug=False)
 
         exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
diff --git a/hysop/methods.py b/hysop/methods.py
index a3bee4ab1b26e61e06f149be7e66d3e099271111..c122b82af912ad97e1a4707ad7f19bb43e002f68 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 PolynomialInterpolation
+from hysop.numerics.interpolation.polynomial import PolynomialInterpolator
 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/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index 2e49a40aa22714d837920052535e9238568377b7..b1a27c2aa5cf17c0f933f751c339f2b6a335030d 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -53,80 +53,6 @@ PolynomialInterpolation = EnumFactory.create('PolynomialInterpolation',
         ])
 
 
-from hysop.constants import SpaceDiscretization
-class SpaceDiscretizationMethod(object):
-    """
-    Operator helper to handle space discretization method.
-    """
-    __default_method = {
-        SpaceDiscretization: 2
-    }
-
-    __available_methods = {
-        SpaceDiscretization: (InstanceOf(int), InstanceOf(SpaceDiscretization)),
-    }
-
-    @classmethod
-    def default_method(cls):
-        dm = super(SpaceDiscretizationMethod, cls).default_method()
-        dm.update(cls.__default_method)
-        return dm
-
-    @classmethod
-    def available_methods(cls):
-        am = super(SpaceDiscretizationMethod, cls).available_methods()
-        am.update(cls.__available_methods)
-        return am
-
-    @debug
-    def handle_method(self,method):
-        super(SpaceDiscretizationMethod, self).handle_method(method)
-        sdm = method.pop(SpaceDiscretization)
-        if isinstance(sdm, int):
-            sd = sdm 
-        elif isinstance(sdm, SpaceDiscretization):
-            assert str(sdm).startswith('FDC')
-            sd = int(str(sdm)[3:])
-        else:
-            raise NotImplementedError(sdm)
-        self.space_discretization_method = sdm
-        self.space_discretization = sd
-
-
-class PolynomialInterpolationMethod(SpaceDiscretizationMethod):
-    """
-    Operator helper to handle polynomial interpolation method.
-    """
-    __default_method = {
-        MultiScaleInterpolation: PolynomialInterpolation.CUBIC_FDC2,
-    }
-
-    __available_methods = {
-        MultiScaleInterpolation: InstanceOf(PolynomialInterpolation),
-    }
-
-    @classmethod
-    def default_method(cls):
-        dm = super(PolynomialInterpolationMethod, cls).default_method()
-        dm.update(cls.__default_method)
-        return dm
-
-    @classmethod
-    def available_methods(cls):
-        am = super(PolynomialInterpolationMethod, cls).available_methods()
-        am.update(cls.__available_methods)
-        return am
-
-    @debug
-    def handle_method(self,method):
-        from hysop.numerics.interpolation.polynomial import PolynomialInterpolator
-        super(PolynomialInterpolationMethod, self).handle_method(method)
-        fd = self.space_discretization
-        pi = method.pop(MultiScaleInterpolation)
-        self.polynomial_interpolation_method = pi
-        self.polynomial_interpolator = PolynomialInterpolator.build_interpolator(pi=pi, fd=fd, dim=self.dim)
-
-
 class PolynomialInterpolator(object):
     @classmethod
     def build_interpolator(cls, pi, dim, fd=None,
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index a73b8e975113cbcd89644d4360d11e87ad72f855..793fb86e4c94720e89d1dbd1513fadaf9e68408e 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -12,6 +12,7 @@ from hysop.tools.numpywrappers import npw
 from hysop.tools.decorators import debug
 from hysop.tools.numerics import find_common_dtype
 from hysop.tools.spectral_utils import SpectralTransformUtils
+from hysop.tools.method_utils import PolynomialInterpolationMethod
 from hysop.fields.continuous_field import Field, ScalarField
 from hysop.parameters.scalar_parameter import ScalarParameter
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
@@ -19,7 +20,6 @@ from hysop.core.graph.node_generator import ComputationalGraphNodeGenerator
 from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
 from hysop.core.memory.memory_request import MemoryRequest
 from hysop.operator.base.spectral_operator import SpectralOperatorBase
-from hysop.numerics.interpolation.polynomial import PolynomialInterpolationMethod
 
 
 class SpatialFilterBase(object):
diff --git a/hysop/tools/method_utils.py b/hysop/tools/method_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..7aa96662d68147142b70c3d8eec04e64355bcd2d
--- /dev/null
+++ b/hysop/tools/method_utils.py
@@ -0,0 +1,78 @@
+
+from hysop.constants import SpaceDiscretization
+from hysop.tools.types import check_instance, to_list, first_not_None, InstanceOf
+from hysop.tools.decorators import debug
+from hysop.numerics.interpolation.polynomial import PolynomialInterpolator, PolynomialInterpolation
+
+class SpaceDiscretizationMethod(object):
+    """
+    Operator helper to handle space discretization method.
+    """
+    __default_method = {
+        SpaceDiscretization: 2
+    }
+
+    __available_methods = {
+        SpaceDiscretization: (InstanceOf(int), InstanceOf(SpaceDiscretization)),
+    }
+
+    @classmethod
+    def default_method(cls):
+        dm = super(SpaceDiscretizationMethod, cls).default_method()
+        dm.update(cls.__default_method)
+        return dm
+
+    @classmethod
+    def available_methods(cls):
+        am = super(SpaceDiscretizationMethod, cls).available_methods()
+        am.update(cls.__available_methods)
+        return am
+
+    @debug
+    def handle_method(self,method):
+        super(SpaceDiscretizationMethod, self).handle_method(method)
+        sdm = method.pop(SpaceDiscretization)
+        if isinstance(sdm, int):
+            sd = sdm 
+        elif isinstance(sdm, SpaceDiscretization):
+            assert str(sdm).startswith('FDC')
+            sd = int(str(sdm)[3:])
+        else:
+            raise NotImplementedError(sdm)
+        self.space_discretization_method = sdm
+        self.space_discretization = sd
+
+
+class PolynomialInterpolationMethod(SpaceDiscretizationMethod):
+    """
+    Operator helper to handle polynomial interpolation method.
+    """
+    __default_method = {
+        PolynomialInterpolator: PolynomialInterpolation.LINEAR,
+    }
+
+    __available_methods = {
+        PolynomialInterpolator: InstanceOf(PolynomialInterpolation),
+    }
+
+    @classmethod
+    def default_method(cls):
+        dm = super(PolynomialInterpolationMethod, cls).default_method()
+        dm.update(cls.__default_method)
+        return dm
+
+    @classmethod
+    def available_methods(cls):
+        am = super(PolynomialInterpolationMethod, cls).available_methods()
+        am.update(cls.__available_methods)
+        return am
+
+    @debug
+    def handle_method(self,method):
+        super(PolynomialInterpolationMethod, self).handle_method(method)
+        fd = self.space_discretization
+        pi = method.pop(PolynomialInterpolator)
+        self.polynomial_interpolation_method = pi
+        self.polynomial_interpolator = PolynomialInterpolator.build_interpolator(pi=pi, fd=fd, dim=self.dim)
+
+