ksdj

HySoP/hysop/operator/discrete/drag_and_lift.py

@@ -321,9 +321,11 @@ class NocaForces(Forces):
             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._buffer, -1,
+                                                         2 * d)
             self.force += self._buffer
-            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer, 1)
+            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer,
+                                                         1, 2 * d + 1)
             self.force += self._buffer
 
         self.force *= self._normalization
@@ -352,9 +354,32 @@ class NocaForces(Forces):
             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)
+                                                         -1, 2 * d)
             self.force += self._buffer * dsurf
-            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer, 1)
+            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer,
+                                                         1, 2 * d + 1)
+            self.force += self._buffer * dsurf
+
+        self.force *= self._normalization
+
+    def _noca_III(self, dt):
+        """
+        Computes local values of the forces using formula 2.10
+        ("Flux Equation") from Noca99
+        @param dt : current time step
+        @return the local (i.e. current mpi process) forces
+        """
+
+        # -- 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_flux(s_x0, self._buffer,
+                                                         -1, 2 * d)
+            self.force += self._buffer * dsurf
+            self._buffer[...] = self._integrate_gamm_imp(s_xn, self._buffer,
+                                                         1, 2 * d + 1)
             self.force += self._buffer * dsurf
 
         self.force *= self._normalization
@@ -439,3 +464,22 @@ class NocaForces(Forces):
             buff[...] = coords[j] * buff[...]
             res[i_n] += npw.real_sum(coords[j] * buff)
             res[j] -= x0 * npw.real_sum(buff)
+
+
+    def _integrate_gamma_flux(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)