diff --git a/HySoP/hysop/domain/subsets/porous.py b/HySoP/hysop/domain/subsets/porous.py
index 25b2454c85b45fe08456eb9dfc89f9474693d761..d6f6aa4c4c4ff48d47ea91f087dc709e294769da 100644
--- a/HySoP/hysop/domain/subsets/porous.py
+++ b/HySoP/hysop/domain/subsets/porous.py
@@ -53,7 +53,7 @@ class Porous(Subset):
                                   radius=in_radius)
             out_set = self._source(parent=self._parent, origin=self.origin,
                                    radius=out_radius)
-            self.ind[topo][i] = Subset.substraction(out_set, in_set, topo)
+            self.ind[topo][i] = Subset.subtract(out_set, in_set, topo)
             out_radius = in_radius
 
         self.ind[topo][-1] = in_set.discretize(topo)
diff --git a/HySoP/hysop/domain/subsets/sphere.py b/HySoP/hysop/domain/subsets/sphere.py
index cc48610b43582e12616e195bdcacafbd9e7ca0a7..290146e0e8bb7e503a0c89afc537e56300b0b616 100644
--- a/HySoP/hysop/domain/subsets/sphere.py
+++ b/HySoP/hysop/domain/subsets/sphere.py
@@ -143,7 +143,7 @@ class Porous(Subset):
                                   radius=in_radius)
             out_set = self._source(parent=self._parent, origin=self.origin,
                                    radius=out_radius)
-            self.ind[topo][i] = Subset.substraction(out_set, in_set, topo)
+            self.ind[topo][i] = Subset.subtract(out_set, in_set, topo)
             out_radius = in_radius
 
         self.ind[topo][-1] = in_set.discretize(topo)
diff --git a/HySoP/hysop/domain/subsets/submesh.py b/HySoP/hysop/domain/subsets/submesh.py
index 885634f521e14830a61332e54f7f4a5c45d0c514..b69e3eaf397fde1a256886e04043952896b81e8f 100644
--- a/HySoP/hysop/domain/subsets/submesh.py
+++ b/HySoP/hysop/domain/subsets/submesh.py
@@ -180,7 +180,7 @@ class SubMesh(object):
         This is only useful when one require the computation
         of the intersection of a regular subset and a sphere-like
         subset.
-        See intersection and substraction methods in Subset class.
+        See intersection and subtract methods in Subset class.
         @param coordinates (like topo.mesh.coords)
         @return : an array of boolean (True if inside)
         """
diff --git a/HySoP/hysop/domain/subsets/subset.py b/HySoP/hysop/domain/subsets/subset.py
index 5225eb7d216b8aee0146d994ca34eaecaab15486..bb01a398d7eeab3efa7eec6b3880e13f91afc0a9 100644
--- a/HySoP/hysop/domain/subsets/subset.py
+++ b/HySoP/hysop/domain/subsets/subset.py
@@ -163,12 +163,12 @@ class Subset(object):
         return c0
 
     @staticmethod
-    def substraction(s1, s2, topo):
+    def subtract(s1, s2, topo):
         """
         @param s1 : first subset
         @param s2 : second subset
         @param topo : a parmepy.mpi.topology.Cartesian
-        @return : substract a subset from another for a
+        @return : subtract a subset from another for a
         given topo and return a tuple of arrays, like self.ind[topo]
         """
         if s1.chi is not None:
@@ -182,12 +182,12 @@ class Subset(object):
         return np.where(np.logical_and(c1, c2))
 
     @staticmethod
-    def substract_list_of_sets(s1, s2, topo):
+    def subtract_list_of_sets(s1, s2, topo):
         """
         @param s1 : a list of subsets
         @param s2 : a second list of subsets
         @param topo : a parmepy.mpi.topology.Cartesian
-        @return : substract a the union of subsets from the union of
+        @return : subtract a the union of subsets from the union of
         other subsets for a
         given topo and return a tuple of arrays, like self.ind[topo]
         """
diff --git a/HySoP/hysop/domain/tests/test_subset.py b/HySoP/hysop/domain/tests/test_subset.py
index 9c5bb949868275bf6f9d361f31350d1df53f40ae..61c63677aada4025eee0ba8d790f31ddb34caaf9 100644
--- a/HySoP/hysop/domain/tests/test_subset.py
+++ b/HySoP/hysop/domain/tests/test_subset.py
@@ -131,7 +131,7 @@ def test_union_2D():
     assert init_obst_list(discr2D, 'multi_obst_2d')
 
 
-def init_substract(discr, fileref):
+def init_subtract(discr, fileref):
     topo, sd, vd = init_multi(discr, fileref)
     dim = len(discr.resolution)
     rd = ldom[0] * 0.1
@@ -142,21 +142,21 @@ def init_substract(discr, fileref):
     s1.discretize(topo)
     box.discretize(topo)
     sd = Field(domain=topo.domain, isVector=False).discretize(topo)
-    indices = Subset.substraction(box, s1, topo)
+    indices = Subset.subtract(box, s1, topo)
     sd[0][indices] = 2.
     iCompute = topo.mesh.iCompute
     return np.allclose(sd[0][iCompute], vd[0][iCompute])
 
 
-def test_substraction_3D():
-    assert init_substract(discr3D, 'multi_obst_subs_3d')
+def test_subtract_3D():
+    assert init_subtract(discr3D, 'multi_obst_subs_3d')
 
 
-def test_substraction_2D():
-    assert init_substract(discr2D, 'multi_obst_subs_2d')
+def test_subtract_2D():
+    assert init_subtract(discr2D, 'multi_obst_subs_2d')
 
 
-def init_substract_lists(discr, fileref):
+def init_subtract_lists(discr, fileref):
     topo, sd, vd = init_multi(discr, fileref)
     dim = len(discr.resolution)
     rd = ldom[0] * 0.1
@@ -180,17 +180,17 @@ def init_substract_lists(discr, fileref):
     for ob in box:
         ob.discretize(topo)
     sd = Field(domain=topo.domain, isVector=False).discretize(topo)
-    indices = Subset.substract_list_of_sets(box, obs, topo)
+    indices = Subset.subtract_list_of_sets(box, obs, topo)
     sd[0][indices] = 2.
     return np.allclose(sd[0][topo.mesh.iCompute], vd[0][topo.mesh.iCompute])
 
 
-def test_substraction_list_3D():
-    assert init_substract_lists(discr3D, 'multi_obst_subslist_3d')
+def test_subtract_list_3D():
+    assert init_subtract_lists(discr3D, 'multi_obst_subslist_3d')
 
 
-def test_substraction_list_2D():
-    assert init_substract_lists(discr2D, 'multi_obst_subslist_2d')
+def test_subtract_list_2D():
+    assert init_subtract_lists(discr2D, 'multi_obst_subslist_2d')
 
 
 # This may be useful to run mpi tests
@@ -203,7 +203,7 @@ if __name__ == "__main__":
     test_hemisphere_2d()
     test_union_2D()
     test_union_3D()
-    test_substraction_3D()
-    test_substraction_2D()
-    test_substraction_list_3D()
-    test_substraction_list_2D()
+    test_subtract_3D()
+    test_subtract_2D()
+    test_subtract_list_3D()
+    test_subtract_list_2D()
diff --git a/HySoP/hysop/numerics/utils.py b/HySoP/hysop/numerics/utils.py
index c07d0c8def6a152f14183333f80abc62953de6c6..30537774a0b6304135660fd354e13bac95079e95 100644
--- a/HySoP/hysop/numerics/utils.py
+++ b/HySoP/hysop/numerics/utils.py
@@ -22,29 +22,9 @@ class Utils(object):
         dim = len(y)
         res = npw.zeros(dim)
         for (i, j) in Utils.gen_indices:
-            work[sl] = x[i] * y[j][sl]
-            work[sl] -= x[j] * y[i][sl]
-            res[current_dir] = npw.real_sum(work[sl])
-            current_dir += 1
-        return res
-
-    @staticmethod
-    def noca_1(x, v, w, sl, work):
-        """
-        @param x : a tuple representing grid coordinates (like coords in mesh)
-        @param y : a discrete field (i.e. list of numpy arrays)
-        @param work : numpy array
-        Sum on a volume (defined by sl) of cross products of
-        x with y at each grid point. Result in work.
-        """
-        current_dir = 0
-        dim = len(w)
-        res = npw.zeros(dim)
-        for (i, j) in Utils.gen_indices:
-            work[sl] = x[i] * w[j][sl]
-            work[sl] -= x[j] * w[i][sl]
-            work[sl] *= v[sl]
-            res[current_dir] = npw.real_sum(work[sl])
+            work[...] = x[i] * y[j][sl]
+            work[...] -= x[j] * y[i][sl]
+            res[current_dir] = npw.real_sum(work[...])
             current_dir += 1
         return res
 
diff --git a/HySoP/hysop/operator/computational.py b/HySoP/hysop/operator/computational.py
index 80b30340ce94ec313c0c7edf0b120769eacf5ec2..8004405ba6cd05c427d8f7acb6017ee0d2f81d5b 100644
--- a/HySoP/hysop/operator/computational.py
+++ b/HySoP/hysop/operator/computational.py
@@ -53,7 +53,7 @@ class Computational(Operator):
         self.method = method
 
         ## The discretization of this operator.
-        self.discreteOperator = None
+        self.discrete_op = None
 
         ## A dictionnary of discreteFields associated with this operator
         ## key = continuous variable \n
@@ -287,8 +287,8 @@ class Computational(Operator):
         """
         Memory cleaning.
         """
-        if self.discreteOperator is not None:
-            self.discreteOperator.finalize()
+        if self.discrete_op is not None:
+            self.discrete_op.finalize()
 
     @debug
     @profile
@@ -300,14 +300,14 @@ class Computational(Operator):
         parmepy.problem.simulation.Simulation for details.
         """
         super(Computational, self).apply(simulation)
-        if self.discreteOperator is not None:
-            self.discreteOperator.apply(simulation)
+        if self.discrete_op is not None:
+            self.discrete_op.apply(simulation)
 
     def printComputeTime(self):
         """ Time monitoring."""
-        if self.discreteOperator is not None:
-            self.discreteOperator.printComputeTime()
-            self.time_info = self.discreteOperator.time_info
+        if self.discrete_op is not None:
+            self.discrete_op.printComputeTime()
+            self.time_info = self.discrete_op.time_info
         else:
             from parmepy.mpi.main_var import main_rank
             shortName = str(self.__class__).rpartition('.')[-1][0:-2]
@@ -320,19 +320,19 @@ class Computational(Operator):
         Update ghost points values, if any.
         """
         assert self._is_discretized
-        self.discreteOperator.update_ghosts()
+        self.discrete_op.update_ghosts()
 
     def __str__(self):
         """
         Common printings for operators
         """
         shortName = str(self.__class__).rpartition('.')[-1][0:-2]
-        if self.discreteOperator is not None:
-            s = str(self.discreteOperator)
+        if self.discrete_op is not None:
+            s = str(self.discrete_op)
         else:
             s = shortName + " operator. Not discretised."
         return s + "\n"
 
     def get_profiling_info(self):
-        if self.discreteOperator is not None:
-            self.profiler += self.discreteOperator.profiler
+        if self.discrete_op is not None:
+            self.profiler += self.discrete_op.profiler
diff --git a/HySoP/hysop/operator/custom.py b/HySoP/hysop/operator/custom.py
index bce5b32f887ba65697dbacbce8a080b56c812239..774e98855d8ad7e8356a8829f49c4ab397b0008e 100644
--- a/HySoP/hysop/operator/custom.py
+++ b/HySoP/hysop/operator/custom.py
@@ -17,9 +17,9 @@ class CustomMonitor(Computational):
     def setup(self, rwork=None, iwork=None):
         if not self._is_uptodate:
             self.discretize()
-            self.discreteOperator = CM(self.function, self.res_shape,
+            self.discrete_op = CM(self.function, self.res_shape,
                                        variables=self.discreteFields.values())
             # Output setup
             self._set_io(self.function.__name__, (1, 1 + self.res_shape))
-            self.discreteOperator.setWriter(self._writer)
+            self.discrete_op.setWriter(self._writer)
             self._is_uptodate = True
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index d0e1267a02e32e8ae65a5ce9f5bdcb68929964c3..3d5c7e7054f9941978bd0d5c4673c0ecf08304f2 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 """
-@file compute_forces.py
-Compute forces
+@file operator/discrete/drag_and_lift.py
+Methods to compute drag and lift forces
 """
 from parmepy.numerics.update_ghosts import UpdateGhosts
 from parmepy.operator.discrete.discrete import DiscreteOperator
@@ -41,7 +41,6 @@ class Forces(DiscreteOperator):
         # deal with obstacles, volume of control ...
         self._indices = self._init_indices(obstacles)
 
-        self._dim = self.domain.dimension
         msg = 'Force computation undefined for domain of dimension 1.'
         assert self._dim > 1, msg
 
@@ -120,6 +119,7 @@ class Forces(DiscreteOperator):
         if not self._on_proc:
             self._buffer[...] = 0.0
             self.force[...] = 0.0
+            self._previous[...] = 0.0
         else:
             self._formula(dt)
         # Reduce results over MPI processes
@@ -131,7 +131,6 @@ class Forces(DiscreteOperator):
             self._writer.buffer[0, 1:] = self.force
             self._writer.write()
 
-
 class MomentumForces(Forces):
     """
     Compute drag and lift using Noca's formula.
@@ -178,7 +177,6 @@ class MomentumForces(Forces):
                                                      component=d)
              for d in xrange(nbc)]
         self.force[...] = 1. / dt * (self._buffer - self._previous)
-
         # Update previous for next time step ...
         self._previous[...] = self._buffer[...]
 
@@ -192,7 +190,7 @@ class NocaForces(Forces):
     """
 
     def __init__(self, velocity, vorticity, nu, volume_of_control,
-                 normalization=1., **kwds):
+                 normalization=1., surfdir=None, **kwds):
         """
         @param vorticity field
         @param nu viscosity
@@ -235,6 +233,14 @@ class NocaForces(Forces):
         self._formula = self._noca
         self._subcomm = self._voc.subcomm[self._topology]
 
+        if surfdir is None:
+            surfdir = [0]
+        # Directions where integrals on surfaces are computed
+        self._surfdir = surfdir
+        # buffer used to compute d/dt velocity
+        # see integrate_gamma_mom
+        self._previous_sum_u = npw.zeros(2 * self._dim)
+
     def _init_indices(self, obstacles):
         """
         Compute a list of indices corresponding to points inside
@@ -254,52 +260,106 @@ class NocaForces(Forces):
                 obs.discretize(self._topology)
             # return indices of points inside the box, excluding points
             # inside the obstacles.
-            return Subset.substract_list_of_sets([self._voc], obstacles,
+            return Subset.subtract_list_of_sets([self._voc], obstacles,
                                                  self._topology)
 
     def _set_work_arrays(self, rwork=None, iwork=None):
-        shape_v = [None, ] * self._dim
+        shape_v = [None, ] * (self._dim + 1)
+        toporef = self.variables[0].topology
         for i in xrange(self._dim):
-            v_ind = self._voc.surf[2 * i].mesh[self._topology].ind4integ
-            shape_v[i] = self.velocity.data[i][v_ind].shape
+            v_ind = self._voc.surf[2 * i].mesh[toporef].ind4integ
+            shape_v[i] = self.velocity.data[0][v_ind].shape
+        v_ind = self._voc.mesh[toporef].ind4integ
+        shape_v[-1] = self.velocity.data[0][v_ind].shape
         # setup for rwork, iwork is useless.
         if rwork is None:
             # ---  Local allocation ---
-            self._rwork = [npw.zeros(shape_v[d]) for d in xrange(self._dim)]
+            self._rwork = [npw.zeros(shape_v[d])
+                           for d in xrange(2 * self._dim)]
+            self._rwork.append(npw.zeros(shape_v[-1]))
         else:
             assert isinstance(rwork, list), 'rwork must be a list.'
             # --- External rwork ---
             self._rwork = rwork
-            assert len(self._rwork) == self._dim
-            for i in xrange(self._dim):
+            assert len(self._rwork) == 2 * self._dim + 1
+            for i in xrange(self._dim + 1):
                 assert self._rwork[i].shape == shape_v[i]
 
-    def _noca(self, dt):
+    def _noca_I(self, dt):
         """
-        Perform integrals on volume and surfaces of the control box
+        Computes local values of the forces using formula 2.1
+        ("Impulse Equation") from Noca99
         @param dt : current time step
-        @return the local (i.e. current mpi process) force
+        @return the local (i.e. current mpi process) forces
         """
         # -- Integration over the volume of control --
         # -1/(N-1) . d/dt int(x ^ w)
-        self._buffer = self._integrateOnBox(self._buffer)
+        mesh = self._voc.mesh[self._topology]
+        coords = mesh.coords4int
+        # slices representing the points inside the subset
+        # (global indices relative to the global mesh of the current topology)
+        sl = mesh.ind4integ
+        if self._dim == 2:
+            wz  = self.vorticity.data[0]
+            self._rwork[-1][...] = coords[YDIR] * wz[sl]
+            self._buffer[0] = npw.real_sum(self._rwork[-1])
+            self._rwork[-1][...] = -coords[XDIR] * wz[sl]
+            self._buffer[1] = npw.real_sum(self._rwork[-1])
+        elif self._dim == 3:
+            self.buffer[...] = Utils.sum_cross_product(coords,
+                                                       self.vorticity.data,
+                                                       sl, self._rwork[-1])
+        self._buffer[...] *= self._dvol
         self.force[...] = -1. / dt * self._coeff * (self._buffer
                                                     - self._previous)
+        # Update previous for next time step ...
+        self._previous[...] = self._buffer[...]
+
+        # -- Integrals on surfaces --
+        # Only on surf. normal to dir in self._surfdir.
+        for d in self._surfdir:
+            s_x0 = self._voc.surf[2 * d]
+            s_xn = self._voc.surf[2 * d + 1]
+            self._buffer[...] = self._integrate_gamm_imp(s_x0,
+                                                         self._buffer, -1)
+            self.force += self._buffer
+            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer, 1)
+            self.force += self._buffer
+
+        self.force *= self._normalization
+
+    def _noca_II(self, dt):
+        """
+        Computes local values of the forces using formula 2.5
+        ("Momentum Equation") from Noca99
+        @param dt : current time step
+        @return the local (i.e. current mpi process) forces
+        """
+        # -- Integration over the volume of control --
+        # -d/dt int(v)
+        nbc = self.velocity.nbComponents
+        self._buffer[...] = \
+            [self._voc.integrate_dfield_on_proc(self.velocity,
+                                                component=d)
+             for d in xrange(nbc)]
+        self.force[...] = -1. / dt * (self._buffer - self._previous)
 
         # Update previous for next time step ...
         self._previous[...] = self._buffer[...]
         # -- Integrals on surfaces --
-        # Only on surf. normal to XDIR
-        s_x0 = self._voc.surf[0]
-        s_xn = self._voc.surf[1]
-        self._buffer = self._integrateOnSurface(s_x0, self._buffer, 1)
-        self.force += self._buffer
-        self._buffer = self._integrateOnSurface(s_xn, self._buffer, -1)
-        self.force += self._buffer
+        # Only on surf. normal to dir in self._surfdir.
+        for d in self._surfdir:
+            s_x0 = self._voc.surf[2 * d]
+            s_xn = self._voc.surf[2 * d + 1]
+            self._buffer[...] = self._integrate_gamm_imp(s_x0, self._buffer,
+                                                         -1)
+            self.force += self._buffer * dsurf
+            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer, 1)
+            self.force += self._buffer * dsurf
 
         self.force *= self._normalization
 
-    def _integrateOnSurface(self, surf, res, normal):
+    def _integrate_gamm_imp(self, surf, res, normal, nsurf):
         res[...] = 0.0
         if not surf.on_proc[self._topology]:
             return res
@@ -330,7 +390,7 @@ class NocaForces(Forces):
         #  1/(dim - 1) * int( (nw)(x X u) - (nu)(x X w))
         wd = self.vorticity.data
         # buffer size == vd[sl].shape, allocated in rwork[i_n]
-        buff = self._rwork[i_n]
+        buff = self._rwork[nsurf]
         # Indices used for cross-product
         j1 = [YDIR, ZDIR, XDIR]
         j2 = [ZDIR, XDIR, YDIR]
@@ -360,21 +420,22 @@ class NocaForces(Forces):
             self._fd_scheme.compute(vd[j], i_n, buff)
             res[j] += normal * self.nu * npw.real_sum(buff)
 
-        res *= dsurf
         return res
 
-    def _integrateOnBox(self, res):
-        """
-        Compute the integral over the volume of control
-        of \f$ x ^ \omega \f$
-        """
-        assert self._dim == 3, 'Not defined for dim < 3'
-        # coords of the points in the volume of control
-        coords = self._voc.coords[self._topology]
-        # slices representing the points inside the subset
-        # (global indices relative to the global mesh of the current topology)
-        sl = self._voc.ind[self._topology]
-        res = Utils.sum_cross_product(coords, self.vorticity.data,
-                                      sl, self._rwork)
-        res *= self._dvol
-        return res
+    def _integrate_gamma_mom(self, surf, res, normal, nsurf):
+        # i_n : normal dir
+        # i_t : other dirs
+        i_n = surf.n_dir
+        i_t = surf.t_dir
+        coords = surf.mesh[self._topology].coords4int
+        x0 = coords[i_n].flat[0]
+        res += self._integrate_gamma_imp(surf, res, normal, nsurf)
+        for j in i_t:
+            buff = self._rwork[nsurf]
+            # buff is supposed to handle velocity at the previous time step
+            buff[...] -= self.velocity.data[j][sl]
+            # compute diff(velocity,t ) in buff
+            buff[...] *= - 1. / dt
+            buff[...] = coords[j] * buff[...]
+            res[i_n] += npw.real_sum(coords[j] * buff)
+            res[j] -= x0 * npw.real_sum(buff)
diff --git a/HySoP/hysop/operator/drag_and_lift.py b/HySoP/hysop/operator/drag_and_lift.py
index 2f97c4b9b8c2d6c1e93ed9aa54bb8c6cf454f835..acddc7ffea79733647a7157afa43383b40067033 100644
--- a/HySoP/hysop/operator/drag_and_lift.py
+++ b/HySoP/hysop/operator/drag_and_lift.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 """
-@file compute_forces.py
-Compute forces
+@file drag_and_lift.py
+Methods to compute drag and lift
 """
 from parmepy.operator.computational import Computational
 from parmepy.operator.continuous import opsetup
@@ -30,14 +30,7 @@ class Forces(Computational):
 
     @abstractmethod
     def get_work_properties(self):
-        if not self._is_discretized:
-            msg = 'The operator must be discretized '
-            msg += 'before any call to this function.'
-            raise RuntimeError(msg)
-        vd = self.discreteFields[self.input[0]]
-        v_ind = vd.topology.mesh.iCompute
-        shape_v = vd[0][v_ind].shape
-        return {'rwork': None, 'iwork': None}
+        pass
 
     def discretize(self):
         super(Forces, self)._standard_discretize(self._min_ghosts)
@@ -53,13 +46,13 @@ class Forces(Computational):
         """
         @return the last computed value for drag force
         """
-        return self.discreteOperator.force[-1]
+        return self.discrete_op.force[0]
 
     def lift(self):
         """
-        @return the last computed value for lift
+        @return the last computed values for lift forces
         """
-        return self.discreteOperator.force[:self.domain.dimension - 1]
+        return self.discrete_op.force[1:]
 
 
 class MomentumForces(Forces):
@@ -90,13 +83,13 @@ class MomentumForces(Forces):
         if not self._is_uptodate:
             from parmepy.operator.discrete.drag_and_lift \
                 import MomentumForces as DMF
-            self.discreteOperator = DMF(
+            self.discrete_op = DMF(
                 velocity=self.discreteFields[self.velocity],
                 obstacles=self.obstacles)
 
             # output setup
             self._set_io('drag_and_lift', (1, 4))
-            self.discreteOperator.setWriter(self._writer)
+            self.discrete_op.setWriter(self._writer)
             self._is_uptodate = True
 
 
@@ -110,7 +103,7 @@ class NocaForces(Forces):
     __metaclass__ = ABCMeta
 
     def __init__(self, velocity, vorticity, nu, volume_of_control=None,
-                 normalization=1., **kwds):
+                 normalization=1., surfdir=None, **kwds):
         """
         @param velocity field
         @param vorticity field
@@ -145,27 +138,35 @@ class NocaForces(Forces):
         elif SpaceDiscretisation is FD_C_4:
             self._min_ghosts = 2
 
+        if surfdir is None:
+            surfdir = [0]
+        # Directions where integrals on surfaces are computed
+        self._surfdir = surfdir
+
         # The volume of control, in which forces are computed
         if volume_of_control is None:
-            xr = self.domain.origin
-            lr = self.domain.length
+            lr = self.domain.length * 0.9
+            xr = self.domain.origin + 0.04 * self.domain.length
         from parmepy.domain.subsets.control_box import ControlBox
         volume_of_control = ControlBox(parent=self.domain,
                                        origin=xr, length=lr)
         self._voc = volume_of_control
 
+
     def get_work_properties(self):
         if not self._is_discretized:
             msg = 'The operator must be discretized '
             msg += 'before any call to this function.'
             raise RuntimeError(msg)
 
-        shape_v = [None, ] * self.domain.dimension
+        shape_v = [None, ] * (self.domain.dimension + 1)
         slist = self._voc.surf
         toporef = self.discreteFields[self.velocity].topology
         for i in xrange(self.domain.dimension):
             v_ind = slist[2 * i].mesh[toporef].ind4integ
             shape_v[i] = self.velocity.data[i][v_ind].shape
+        v_ind = self._voc.mesh[toporef].ind4integ
+        shape_v[-1] = self.velocity.data[0][v_ind].shape
         # setup for rwork, iwork is useless.
         return {'rwork': shape_v, 'iwork': None}
 
@@ -176,15 +177,16 @@ class NocaForces(Forces):
                 import NocaForces as DNF
             topo = self.discreteFields[self.velocity].topology
             self._voc.discretize(topo)
-            self.discreteOperator = DNF(
+            self.discrete_op = DNF(
                 velocity=self.discreteFields[self.velocity],
                 vorticity=self.discreteFields[self.vorticity],
                 nu=self.nu,
                 volume_of_control=self._voc,
                 normalization=self.normalization,
-                obstacles=self.obstacles)
+                obstacles=self.obstacles,
+                surfdir=self._surfdir)
 
             # output setup
             self._set_io('drag_and_lift', (1, 4))
-            self.discreteOperator.setWriter(self._writer)
+            self.discrete_op.setWriter(self._writer)
             self._is_uptodate = True
diff --git a/HySoP/hysop/operator/penalization_and_curl.py b/HySoP/hysop/operator/penalization_and_curl.py
index cc3e6b9a5ce50c36d8e7572c61b6fe153ce6306c..537b7b31ec61b3a825b785e22f95fb82e83fbe61 100644
--- a/HySoP/hysop/operator/penalization_and_curl.py
+++ b/HySoP/hysop/operator/penalization_and_curl.py
@@ -78,9 +78,9 @@ class PenalizationAndCurl(Penalization):
     @opsetup
     def setup(self, rwork=None, iwork=None):
         self._curl.setup(rwork, iwork)
-        self.discreteOperator = PAC(
+        self.discrete_op = PAC(
             vorticity=self.discreteFields[self.vorticity],
-            curl=self._curl.discreteOperator,
+            curl=self._curl.discrete_op,
             obstacles=self.obstacles,
             coeff=self.coeff, rwork=rwork, iwork=iwork)
         self._is_uptodate = True
diff --git a/HySoP/hysop/operator/tests/test_adaptive_time_step_flymake/p1/dt_adapt b/HySoP/hysop/operator/tests/test_adaptive_time_step_flymake/p1/dt_adapt
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/HySoP/hysop/operator/tests/test_drag_and_lift.py b/HySoP/hysop/operator/tests/test_drag_and_lift.py
index c7645fa2e7570767d56f55f3b7709f7132726acd..89c315480004e2f45bf161ebaac25161542dc8a3 100644
--- a/HySoP/hysop/operator/tests/test_drag_and_lift.py
+++ b/HySoP/hysop/operator/tests/test_drag_and_lift.py
@@ -99,7 +99,7 @@ def test_momentum_forces_3D():
     """
     topo, wref, vref = init_vort(discr3D, 'drag_3d_sphere')
     # Obstacles
-    rd = ldom[0] * 0.3
+    rd = ldom[0] * 0.2
     # Fields to penalize
     vorti = Field(domain=topo.domain, formula=v3d, name='Vorti', isVector=True)
     velo = Field(domain=topo.domain, formula=v3d, name='Velo', isVector=True)
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D_.xmf b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D_.xmf
deleted file mode 100644
index 5a79944cd588471e28e65fce8ddafab0611a109f..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D_.xmf
+++ /dev/null
@@ -1,59 +0,0 @@
-<?xml version="1.0" ?>
-<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd">
-<Xdmf Version="2.0">
- <Domain>
-  <Grid Name="CellTime" GridType="Collection" CollectionType="Temporal">
-   <Grid Name="Iteration 000" GridType="Uniform">
-    <Time Value="0.333333333333" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="64  64  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     -1.0  -1.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.03125  0.03125  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="ScalRef3D_1_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="64  64  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      ScalRef3D__00000.h5:/ScalRef3D_1_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 001" GridType="Uniform">
-    <Time Value="0.666666666667" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="64  64  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     -1.0  -1.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.03125  0.03125  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="ScalRef3D_1_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="64  64  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      ScalRef3D__00001.h5:/ScalRef3D_1_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 002" GridType="Uniform">
-    <Time Value="1.0" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="64  64  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     -1.0  -1.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.03125  0.03125  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="ScalRef3D_1_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="64  64  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      ScalRef3D__00002.h5:/ScalRef3D_1_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-  </Grid>
- </Domain>
-</Xdmf>
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00000.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00000.h5
deleted file mode 100644
index 70243a0dfc538b43f9c0e77b388e7c39d70b7ef4..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00000.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00001.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00001.h5
deleted file mode 100644
index 70243a0dfc538b43f9c0e77b388e7c39d70b7ef4..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00001.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00002.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00002.h5
deleted file mode 100644
index 70243a0dfc538b43f9c0e77b388e7c39d70b7ef4..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/ScalRef3D__00002.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal.xmf b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal.xmf
deleted file mode 100644
index 167c44c016c1bae9db38c54b879c7400eaa3d5f1..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal.xmf
+++ /dev/null
@@ -1,42 +0,0 @@
-<?xml version="1.0" ?>
-<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd">
-<Xdmf Version="2.0">
- <Domain>
-  <Grid Name="CellTime" GridType="Collection" CollectionType="Temporal">
-   <Grid Name="Iteration 000" GridType="Uniform">
-    <Time Value="0.0" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="64  64  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     -1.0  -1.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.03125  0.03125  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="ScalRef3D_0_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="64  64  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      testIO_scal_00000.h5:/ScalRef3D_0_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 001" GridType="Uniform">
-    <Time Value="0.3" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="64  64  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     -1.0  -1.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.03125  0.03125  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="ScalRef3D_0_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="64  64  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      testIO_scal_00001.h5:/ScalRef3D_0_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-  </Grid>
- </Domain>
-</Xdmf>
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal_00000.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal_00000.h5
deleted file mode 100644
index 286edb5df96691b9883828f22b53c311d1ba3f24..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal_00000.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal_00001.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal_00001.h5
deleted file mode 100644
index 286edb5df96691b9883828f22b53c311d1ba3f24..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/testIO_scal_00001.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti_.xmf b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti_.xmf
deleted file mode 100644
index 3f1c020b1356959c80ce6649726c0aeb9350ab51..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti_.xmf
+++ /dev/null
@@ -1,134 +0,0 @@
-<?xml version="1.0" ?>
-<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd">
-<Xdmf Version="2.0">
- <Domain>
-  <Grid Name="CellTime" GridType="Collection" CollectionType="Temporal">
-   <Grid Name="Iteration 000" GridType="Uniform">
-    <Time Value="0.333333333333" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="72  129  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     3.8999999999999999  2.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.10833333333333334  0.031007751937984496  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="velo_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00000.h5:/velo_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00000.h5:/velo_2_X
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00000.h5:/velo_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00000.h5:/vorti_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00000.h5:/vorti_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00000.h5:/vorti_2_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 001" GridType="Uniform">
-    <Time Value="0.666666666667" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="72  129  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     3.8999999999999999  2.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.10833333333333334  0.031007751937984496  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="velo_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00001.h5:/velo_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00001.h5:/velo_2_X
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00001.h5:/velo_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00001.h5:/vorti_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00001.h5:/vorti_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00001.h5:/vorti_2_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 002" GridType="Uniform">
-    <Time Value="1.0" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="72  129  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     3.8999999999999999  2.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.10833333333333334  0.031007751937984496  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="velo_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00002.h5:/velo_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00002.h5:/velo_2_X
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00002.h5:/velo_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00002.h5:/vorti_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00002.h5:/vorti_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      velo_vorti__00002.h5:/vorti_2_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-  </Grid>
- </Domain>
-</Xdmf>
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00000.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00000.h5
deleted file mode 100644
index 4840283fbfb035c18e3430a02e61d10f7b96c06f..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00000.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00001.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00001.h5
deleted file mode 100644
index 51b2234556750f019ae940f5f13dc60a7ecef63c..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00001.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00002.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00002.h5
deleted file mode 100644
index ab2b0b2ec34470101e75459841fd9b7df1808254..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/velo_vorti__00002.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo_.xmf b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo_.xmf
deleted file mode 100644
index dfc8f1ed658e69d3bd66fc3310b759f8192f1541..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo_.xmf
+++ /dev/null
@@ -1,134 +0,0 @@
-<?xml version="1.0" ?>
-<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd">
-<Xdmf Version="2.0">
- <Domain>
-  <Grid Name="CellTime" GridType="Collection" CollectionType="Temporal">
-   <Grid Name="Iteration 000" GridType="Uniform">
-    <Time Value="0.333333333333" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="72  129  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     3.8999999999999999  2.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.10833333333333334  0.031007751937984496  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="velo_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00000.h5:/velo_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00000.h5:/velo_2_X
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00000.h5:/velo_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00000.h5:/vorti_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00000.h5:/vorti_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00000.h5:/vorti_2_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 001" GridType="Uniform">
-    <Time Value="0.666666666667" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="72  129  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     3.8999999999999999  2.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.10833333333333334  0.031007751937984496  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="velo_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00001.h5:/velo_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00001.h5:/velo_2_X
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00001.h5:/velo_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00001.h5:/vorti_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00001.h5:/vorti_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00001.h5:/vorti_2_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-   <Grid Name="Iteration 002" GridType="Uniform">
-    <Time Value="1.0" />
-    <Topology TopologyType="3DCORECTMesh" NumberOfElements="72  129  64 "/>
-    <Geometry GeometryType="ORIGIN_DXDYDZ">
-     <DataItem Dimensions="3 " NumberType="Float" Precision="4" Format="XML">
-     3.8999999999999999  2.0  -1.0
-     </DataItem>
-     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
-     0.10833333333333334  0.031007751937984496  0.03125
-     </DataItem>
-    </Geometry>
-    <Attribute Name="velo_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00002.h5:/velo_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00002.h5:/velo_2_X
-     </DataItem>
-    </Attribute>
-    <Attribute Name="velo_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00002.h5:/velo_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Z" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00002.h5:/vorti_2_Z
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_Y" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00002.h5:/vorti_2_Y
-     </DataItem>
-    </Attribute>
-    <Attribute Name="vorti_2_X" AttributeType="Scalar" Center="Node">
-     <DataItem Dimensions="72  129  64 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
-      vorti_velo__00002.h5:/vorti_2_X
-     </DataItem>
-    </Attribute>
-   </Grid>
-  </Grid>
- </Domain>
-</Xdmf>
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00000.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00000.h5
deleted file mode 100644
index d541a445fbc9c6f039fc56f5c63e0978eb68111d..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00000.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00001.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00001.h5
deleted file mode 100644
index d541a445fbc9c6f039fc56f5c63e0978eb68111d..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00001.h5 and /dev/null differ
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00002.h5 b/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00002.h5
deleted file mode 100644
index 22a8bdf14614e8a331786fa00acbb75f185d86be..0000000000000000000000000000000000000000
Binary files a/HySoP/hysop/operator/tests/test_hdf5_io_flymake/p1/vorti_velo__00002.h5 and /dev/null differ