diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 561362e48d14f939661027ab472d0eda3e877393..800fdd02468aa3a881b42bdf83c1d5921291faf3 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -134,6 +134,8 @@ class DragAndLift(Monitoring):
         dt = simulation.timeStep
         ite = simulation.currentIteration
 
+        self.force = npw.zeros(self._dim)
+
         # Synchro of ghost points is required for fd schemes
         self._synchronize(self.vd.data + self.wd.data)
 
@@ -208,7 +210,7 @@ class DragAndLift(Monitoring):
             res[j] -= self._coeff * normal[i1] * coords[i1] * np.sum(buff)
 
         # Last part
-        # Update fd schemes so to compute laplacian and other derivatives
+        # Update fd schemes in order to compute laplacian and other derivatives
         # only on the surface (i.e. for list of indices in sl)
         self._laplacian.fd_scheme.computeIndices(sl)
         for j in i2:
@@ -219,10 +221,12 @@ class DragAndLift(Monitoring):
                 np.sum(self._work[sl])
         self._fd_scheme.computeIndices(sl)
         self._fd_scheme.compute(vdata[i1], i1, self._work)
-        res[i1] += normal[i1] * self.nu * np.sum(self._work[sl])
+        res[i1] += 2.0 * normal[i1] * self.nu * np.sum(self._work[sl])
         for j in i2:
             self._fd_scheme.compute(vdata[i1], j, self._work)
             res[j] += normal[i1] * self.nu * np.sum(self._work[sl])
+            self._fd_scheme.compute(vdata[j], i1, self._work)
+            res[j] += normal[i1] * self.nu * np.sum(self._work[sl])
 
         res *= dsurf
         return res