diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index b639d037d376f4965af8db22736552d7ca7c5109..b243b260e625b62926f7a536543ac434a61a5267 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -31,13 +31,14 @@ CartesianGrid = grid.CartesianGrid
 
 
 ## Fields
-#import fields.discrete
-#import fields.continuous
-#import fields.analytical
-#ScalarField = fields.discrete.ScalarField
-#VectorField = fields.discrete.VectorField
-#ContinuousField = fields.continuous.ContinuousField
-#AnalyticalField = fields.analytical.AnalyticalField
+import fields.discrete
+import fields.continuous
+import fields.analytical
+ScalarField = fields.discrete.ScalarField
+VectorField = fields.discrete.VectorField
+ContinuousField = fields.continuous.ContinuousField
+AnalyticalField = fields.analytical.AnalyticalField
+
 
 ## Obstacle
 ## import obstacle.obstacle
@@ -46,25 +47,25 @@ CartesianGrid = grid.CartesianGrid
 ## Obstacle_d = obstacle.obstacle_d.Obstacle_d
 
 ## Operators
-### import operator.advec_scales
+## import operator.advec_scales
 ### import operator.diffusion
 ## import operator.penalization
 ## ### import operator.poisson
-## import operator.stretching
+import operator.stretching
 ## import operator.transport
 ## import operator.velocity
 ## # Advec_scales = operator.advec_scales.Advec_scales
 ## # Diffusion = operator.diffusion.Diffusion
 ## Penalization = operator.penalization.Penalization
 ## # Poisson = operator.poisson.Poisson
-## Stretching = operator.stretching.Stretching
+Stretching = operator.stretching.Stretching
 ## Transport = operator.transport.Transport
 ## Velocity = operator.velocity.Velocity
 
 
 ## ## Problem
-## import problem.problem
-## Problem = problem.problem.Problem
+import problem.problem
+Problem = problem.problem.Problem
 
 
 ## ## Forces
@@ -72,9 +73,9 @@ CartesianGrid = grid.CartesianGrid
 ## Compute_forces = physics.compute_forces.Compute_forces
 
 ## ## Solver
-## # import particular_solvers.basic
+import particular_solvers.basic
 ## #import particular_solvers.gpu
-## # ParticleSolver = particular_solvers.basic.ParticleSolver
+ParticleSolver = particular_solvers.basic.ParticleSolver
 ## #GPUParticleSolver = particular_solvers.gpu.GPUParticleSolver
 
 ## ## Tools
diff --git a/HySoP/hysop/integrator/__init__.py b/HySoP/hysop/integrator/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..4a661fed99196f000dc6530c4bf4f53868907bc4
--- /dev/null
+++ b/HySoP/hysop/integrator/__init__.py
@@ -0,0 +1,15 @@
+"""
+@package parmepy.particular_solvers.integrator
+
+Everything concerning time integration.
+
+"""
+import runge_kutta2
+import runge_kutta3
+import runge_kutta4
+RK2 = runge_kutta2.RK2
+RK3 = runge_kutta3.RK3
+RK4 = runge_kutta4.RK4
+
+import euler
+euler = euler.Euler
diff --git a/HySoP/hysop/particular_solvers/integrator/euler.py b/HySoP/hysop/integrator/euler.py
similarity index 73%
rename from HySoP/hysop/particular_solvers/integrator/euler.py
rename to HySoP/hysop/integrator/euler.py
index 819c94712b6720acafad7c229a7bd90191f74d54..4e58120343b3b186580401e4214e7d38f3a9ef36 100755
--- a/HySoP/hysop/particular_solvers/integrator/euler.py
+++ b/HySoP/hysop/integrator/euler.py
@@ -1,18 +1,18 @@
 # -*- coding: utf-8 -*-
 """
-@package Utils
+@package parmepy.ode_solver.Euler
+
 Forward Euler method.
 """
-from .integrator import ODESolver
+from ode_solver import ODESolver
 import numpy as np
-
-
-#providedClass = "Euler"
+from parmepy.constants import PARMES_REAL, ORDER
 
 
 class Euler(ODESolver):
     """
-    ODESolver implementation for solving an equation system with forward Euler method.
+    ODESolver implementation for solving an equation system
+    with forward Euler method.
     y'(t) = f(t,y)
 
     y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
@@ -38,10 +38,9 @@ class Euler(ODESolver):
         @param t : current time.
         @param dt : time step.
 #        """
-        result = y + dt * np.asarray(f(t,y))
+        result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
         return result
 
-
     def __str__(self):
         """ToString method"""
         return "Euler Method"
@@ -49,5 +48,5 @@ class Euler(ODESolver):
 
 if __name__ == "__main__":
     print __doc__
-    print "- Provided class : " + providedClass
-    print eval(providedClass).__doc__
+    print "- Provided class : Euler"
+    print Euler.__doc__
diff --git a/HySoP/hysop/particular_solvers/integrator/integrator.py b/HySoP/hysop/integrator/ode_solver.py
similarity index 62%
rename from HySoP/hysop/particular_solvers/integrator/integrator.py
rename to HySoP/hysop/integrator/ode_solver.py
index f418635a0904bc31df6df26b7c5e3aaced35d7e0..33125ba11320c8b9262ba16cca4ceebe3c796999 100644
--- a/HySoP/hysop/particular_solvers/integrator/integrator.py
+++ b/HySoP/hysop/integrator/ode_solver.py
@@ -1,18 +1,25 @@
 # -*- coding: utf-8 -*-
 """
-@package Utils
-ODESolver interface.
+@package parmepy.integrator.ode_solver
+
+Classes to describe (time) integrators.
 """
-from ...constants import *
+from abc import ABCMeta, abstractmethod
 import numpy as np
 
+
 class ODESolver:
     """
-    ODESolver interface for solving an equation system :
+    Abstract description for ODE solvers.
+    Solve the system :
     y'(t) = f(t,y)
 
     By default, f = 0.
     """
+
+    __metaclass__ = ABCMeta
+
+    @abstractmethod
     def __init__(self, f=None):
         """
         Constructor.
@@ -20,11 +27,13 @@ class ODESolver:
         @param f function f(t,y) : Right hand side of the equation to solve.
         """
         ## RHS.
-        if f != None:
+        if f is not None:
             self.f = f
         else:
-            self.f = lambda t, y: np.vectorize(f)(t,y)# [0. for x in y]
+            #[0. for x in y]
+            self.f = lambda t, y: np.vectorize(f)(t, y)
 
+    @abstractmethod
     def integrate(self, y, t, dt):
         """
         Abstract method, apply operaton on a variable.
@@ -35,9 +44,8 @@ class ODESolver:
         @param dt : time step.
         @return y at t+dt.
         """
-        raise NotImplementedError("Need to override method in a subclass of " + providedClass)
 
 if __name__ == "__main__":
     print __doc__
-    print "- Provided class : ODESolver"
+    print "- Provided class : ODESolver (abstract)"
     print ODESolver.__doc__
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py b/HySoP/hysop/integrator/runge_kutta2.py
similarity index 66%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
rename to HySoP/hysop/integrator/runge_kutta2.py
index 6bff109502efb91c1a665bc4be1e7a326c1fc8fa..f6612e7a139800a666ef0f2ac061b277017a62c3 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
+++ b/HySoP/hysop/integrator/runge_kutta2.py
@@ -1,9 +1,10 @@
 # -*- coding: utf-8 -*-
 """
-@package Utils
+@package ode_solver.runge_kutta2
+
 RK2 method interface.
 """
-from .integrator import ODESolver
+from ode_solver import ODESolver
 import numpy as np
 
 
@@ -24,21 +25,22 @@ class RK2(ODESolver):
         """
         ODESolver.__init__(self, f)
 
-    def integrate(self, f, t , dt , y):
+    def integrate(self, f, t, dt, y):
 
         """
         Integration step for RK2 Method.
-         up = f[t, Y(t,x,y,z)].
+         up = f[t, Y(t,space)].
 
-        @param y : function of position at time t.
-        @param f : function of function of position at time t.
+        @param f : right-hand side.
         @param t : current time.
         @param dt : time step.
+        @param y : function Y(t,space)
         """
-#K1 = dt *f[t,y(t,x,y,z)]
-        K1 = dt * np.asarray(f(t,y))
-#result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-        result = y + dt *np.asarray(f(t+dt/2.,y + K1/2.))
+
+        # K1 = dt *f[t,y(t,x,y,z)]
+        K1 = dt * np.asarray(f(t, y))
+        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+        result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.))
         return result
 
     def __str__(self):
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py b/HySoP/hysop/integrator/runge_kutta3.py
similarity index 57%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
rename to HySoP/hysop/integrator/runge_kutta3.py
index dc1d95beab3c67826455acd1df6eae5348c12224..bf51ab41f4124639a78b32d2f18e6d8cf603f82e 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
+++ b/HySoP/hysop/integrator/runge_kutta3.py
@@ -1,9 +1,11 @@
 # -*- coding: utf-8 -*-
 """
-@package Utils
+@package ode_solver.runge_kutta3
+
 RK3 method interface.
 """
-from .integrator import ODESolver
+
+from ode_solver import ODESolver
 import numpy as np
 
 
@@ -25,25 +27,23 @@ class RK3(ODESolver):
         """
         ODESolver.__init__(self, f)
 
-
-    def integrate(self, f, t , dt , y):
-
+    def integrate(self, f, t, dt, y):
         """
-        Integration step for RK3 Method.
-         up = f[t, Y(t,x,y,z)].
+        Integration step for RK2 Method.
+         up = f[t, Y(t,space)].
 
-        @param y : function of position at time t.
-        @param f : function of function of position at time t.
+        @param f : right-hand side.
         @param t : current time.
         @param dt : time step.
+        @param y : function Y(t,space)
         """
-#K1 = dt *f[t,y(t,x,y,z)]
-        K1 = dt * np.asarray(f(t,y))
-#result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-        K2 = dt *np.asarray(f(t+dt/3.,y + K1[:,:,:,:]/3.))
-        K3 = dt *np.asarray(f(t+dt*2./3.,y + K2[:,:,:,:]*2./3.))
-#result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
-        result = y + 1./4.*(K1+3.*K3)
+        # K1 = dt *f[t,y(t,x,y,z)]
+        K1 = dt * np.asarray(f(t, y))
+        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+        K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.))
+        K3 = dt * np.asarray(f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.))
+        # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
+        result = y + 1. / 4. * (K1 + 3. * K3)
         return result
 
     def __str__(self):
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py b/HySoP/hysop/integrator/runge_kutta4.py
similarity index 53%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
rename to HySoP/hysop/integrator/runge_kutta4.py
index 26bdba75037540614967e3da5b7ba0085143d90f..0debd1c1a44b1fadca7cb4b4bb7917e8e99d8f91 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
+++ b/HySoP/hysop/integrator/runge_kutta4.py
@@ -1,9 +1,11 @@
 # -*- coding: utf-8 -*-
 """
-@package Utils
+@package ode_solver.runge_kutta2
+
 RK4 method interface.
 """
-from .integrator import ODESolver
+
+from ode_solver import ODESolver
 import numpy as np
 
 
@@ -26,23 +28,25 @@ class RK4(ODESolver):
         ODESolver.__init__(self, f)
 
 
-    def integrate(self, f, t , dt , y):
+    def integrate(self, f, t, dt, y):
         """
-        Integration step for RK4 Method.
-         up = f[t, Y(t,x,y,z)].
+        Integration step for RK2 Method.
+         up = f[t, Y(t,space)].
 
-        @param y : function of position at time t.
-        @param f : function of function of position at time t.
+        @param f : right-hand side.
         @param t : current time.
         @param dt : time step.
+        @param y : function Y(t,space)
         """
-#K1 = dt *f[t,y(t,x,y,z)]
-        K1 = dt * np.asarray(f(t,y))
-        K2 = dt *np.asarray(f(t+dt/2.,y + K1[:,:,:,:]/2.))
-        K3 = dt *np.asarray(f(t+dt/2.,y + K2[:,:,:,:]/2.))
-        K4 = dt *np.asarray(f(t+dt,y + K3[:,:,:,:]))
-#result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-        result = y + 1./6. * (K1[:,:,:,:] + 2.*K2[:,:,:,:] + 2.*K3[:,:,:,:]+ K4[:,:,:,:])
+
+        # K1 = dt *f[t, y(t, x, y, z)]
+        K1 = dt * np.asarray(f(t, y))
+        K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.))
+        K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.))
+        K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]))
+        # result = y(t, x, y, z)]  +  dt *f[t + dt/2 ,  y(t, x, y, z) + k1/2.]
+        result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
+                                2. * K3[:, :, :, :] + K4[:, :, :, :])
         return result
 
     def __str__(self):
diff --git a/HySoP/hysop/mpi/tests/test_topology.py b/HySoP/hysop/mpi/tests/test_topology.py
index 46ffd526251e2a314ffed564689e51eb42c367a3..b805313e1be7cd9c346c0ef38b5f296518b06106 100644
--- a/HySoP/hysop/mpi/tests/test_topology.py
+++ b/HySoP/hysop/mpi/tests/test_topology.py
@@ -1,6 +1,14 @@
 import parmepy as pp
 
 
+## @pytest.fixture
+## def topology():
+##     import parmepy.mpi.topology as pptopo
+##     import parmepy as pp
+##     #dom = pp.Box()
+##     return resolTopo
+
+
 def test_create_topology1D():
     dom = pp.Box()
     resolTopo = [33, 33, 17]
@@ -9,7 +17,6 @@ def test_create_topology1D():
     assert topo.domain == dom
     assert topo.dim == 1
     assert topo.size == pp.mpi.main_size
-
     print topo
 
 
@@ -44,3 +51,12 @@ def test_create_topology_with_dims():
     assert topo.size == pp.mpi.main_size
     testList = topo.dims == topoDims
     assert testList.all()
+
+
+if __name__ == "__main__":
+
+    #test_create_topology1D()
+    #pp.mpi.main_comm.barrier()
+
+    test_create_topology2D()
+    pp.mpi.main_comm.barrier()
diff --git a/HySoP/hysop/mpi/topology.py b/HySoP/hysop/mpi/topology.py
index b51e6fe0c951f49497106b923ac737edce406af4..27df9147b5547e35733323ea62925926a08c45c7 100644
--- a/HySoP/hysop/mpi/topology.py
+++ b/HySoP/hysop/mpi/topology.py
@@ -136,7 +136,7 @@ class Cartesian(object):
     """
 
     ## Counter of topology.Cartesian instances to set idTopo.
-    _ids = count(0)
+    __ids = count(0)
 
     def __init__(self, domain, dim, globalMeshResolution,
                  comm=main_comm, periods=None, ghosts=None, cutdir=None):
@@ -144,7 +144,7 @@ class Cartesian(object):
         ## Associated domain
         self.domain = domain
         ## An id for the topology
-        self.idTopo = self._ids.next()
+        self.idTopo = self.__ids.next()
         ## Topology dimension
         self.dim = dim
         ## Communicator used to build the topology
@@ -222,7 +222,7 @@ class Cartesian(object):
         ##  Global Mesh Resolution
         self.globalMeshResolution = np.asarray(globalMeshResolution,
                                                dtype=PARMES_INTEGER)
-        
+
         # set cutdir if required (i.e. not set by withResolution constructor).
         if cutdir is None:  # not hasattr(self, '__cutdir'):
             self.__cutdir = np.zeros((domain.dimension), dtype=bool)
@@ -231,7 +231,7 @@ class Cartesian(object):
                 self.__cutdir[:] = self.__cutdir[::-1]
         else:
             self.__cutdir = np.asarray(cutdir, dtype=bool)
-            
+
         # Number of "computed" points (i.e. excluding ghosts/boundaries).
         pts_noghost = np.zeros((self.dim), dtype=PARMES_INTEGER)
         # Warning FP : we should test boundary conditions type here
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index 7f47c6d99730f37d2cffc703585d40891324a191..c18d9bc62331e5e0f2092fa27d46be57aaef6f55 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -15,7 +15,7 @@ from parmepy.operator.poisson import Poisson
 from parmepy.operator.splitting import Splitting
 from parmepy.tools.timer import Timer
 from parmepy.tools.printer import Printer
-
+from parmepy.integrator import RK3
 
 class ParticleSolver(Solver):
     """
diff --git a/HySoP/hysop/particular_solvers/integrator/__init__.py b/HySoP/hysop/particular_solvers/integrator/__init__.py
deleted file mode 100644
index 48c83c107ba4b8011ac4d47b2f8275f204530750..0000000000000000000000000000000000000000
--- a/HySoP/hysop/particular_solvers/integrator/__init__.py
+++ /dev/null
@@ -1,6 +0,0 @@
-"""
-@package parmepy.particular_solvers.integrator
-
-Everything concerning time integration.
-
-"""
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index fce59388f3b7326dee9ce3c4511c1e46ce3a6681..b6109a1d9f5de0b06aa995f52c8c3151f25daeaf 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -3,8 +3,8 @@
 
 Problem representation.
 """
-from ..particular_solvers.basic import ParticleSolver
-from ..particular_solvers.gpu import GPUParticleSolver
+from parmepy.particular_solvers.basic import ParticleSolver
+from parmepy.particular_solvers.gpu import GPUParticleSolver
 
 
 class Problem():
diff --git a/HySoP/hysop/particular_solvers/integrator/forward_euler.py b/HySoP/hysop/unusedOrObsolet/forward_euler.py
similarity index 100%
rename from HySoP/hysop/particular_solvers/integrator/forward_euler.py
rename to HySoP/hysop/unusedOrObsolet/forward_euler.py
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta.py b/HySoP/hysop/unusedOrObsolet/runge_kutta.py
similarity index 100%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta.py
rename to HySoP/hysop/unusedOrObsolet/runge_kutta.py
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py b/HySoP/hysop/unusedOrObsolet/runge_kutta2stretching.py
similarity index 100%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
rename to HySoP/hysop/unusedOrObsolet/runge_kutta2stretching.py
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta3stretching.py b/HySoP/hysop/unusedOrObsolet/runge_kutta3stretching.py
similarity index 100%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta3stretching.py
rename to HySoP/hysop/unusedOrObsolet/runge_kutta3stretching.py
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py b/HySoP/hysop/unusedOrObsolet/runge_kutta4stretching.py
similarity index 100%
rename from HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
rename to HySoP/hysop/unusedOrObsolet/runge_kutta4stretching.py