re-initialize force to zero at the begin of apply --> fix most of the...

re-initialize force to zero at the begin of apply --> fix most of the parallelism problems, but there are still some discrepencies in the solution depending on the number of proc used. These discrepencies are due to the last term of the surface integrals, involving the deformation tensor. Moreover, the last contribution of these surface integrals has been corrected (some terms where missing) 

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