diff --git a/HySoP/hysop/numerics/finite_differences.py b/HySoP/hysop/numerics/finite_differences.py
index bc05d712a6ce3ab699c29d48e5ed014bb1a53007..d90cbec11caeea71c7914499000ebd77622dcf28 100644
--- a/HySoP/hysop/numerics/finite_differences.py
+++ b/HySoP/hysop/numerics/finite_differences.py
@@ -82,6 +82,16 @@ class FiniteDifference(object):
mesh on which finite-differences will be applied.
"""
+ def computeIndices_reduced(self, indices, reduced_shape):
+ """
+ @param indices : a list of slices (see for example
+ parmepy.mpi.mesh.Mesh iCompute) that represent the local
+ mesh on which finite-differences will be applied.
+ """
+ self.computeIndices(indices)
+ self.indices = [slice(0, reduced_shape[d])
+ for d in xrange(len(reduced_shape))]
+
@abstractmethod
def compute(self, tab, cdir, result):
"""
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index 2cd8b4afc892c3a64ea4da2a5ee173451f0d412e..d0e1267a02e32e8ae65a5ce9f5bdcb68929964c3 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -33,6 +33,11 @@ class Forces(DiscreteOperator):
control
"""
super(Forces, self).__init__(**kwds)
+ # topology common to all variables
+ self.input = self.variables
+ self._topology = self.input[0].topology
+ assert (self._topology.ghosts() >= 1).all()
+
# deal with obstacles, volume of control ...
self._indices = self._init_indices(obstacles)
@@ -46,10 +51,6 @@ class Forces(DiscreteOperator):
## The force we want to compute (lift(s) and drag)
self.force = npw.zeros(self._dim)
- # topology common to all variables
- self._topology = self.input[0].topology
- assert (self._topology.ghosts() >= 1).all()
-
# list of np arrays to be synchronized
self._datalist = []
for v in self.input:
@@ -65,7 +66,6 @@ class Forces(DiscreteOperator):
# true if the operator needs to work on the current process.
# Updated in derived class.
self._on_proc = True
- self.input = self.variables
# Set how reduction will be performed
# Default = reduction over all process.
# \todo : add param to choose this option
@@ -88,19 +88,7 @@ class Forces(DiscreteOperator):
pass
def _set_work_arrays(self, rwork=None, iwork=None):
-
- v_ind = self.input[0].topology.mesh.iCompute
- shape_v = self.input[0].data[0][v_ind].shape
- # setup for rwork, iwork is useless.
- if rwork is None:
- # --- Local allocation ---
- self._rwork = [npw.zeros(shape_v)]
- else:
- assert isinstance(rwork, list), 'rwork must be a list.'
- # --- External rwork ---
- self._rwork = rwork
- assert len(self._rwork) == 1
- assert self._rwork[0].shape == shape_v
+ pass
def _mpi_allsum(self):
"""
@@ -270,19 +258,21 @@ class NocaForces(Forces):
self._topology)
def _set_work_arrays(self, rwork=None, iwork=None):
-
- v_ind = self.velocity.topology.mesh.iCompute
- shape_v = self.velocity.data[0][v_ind].shape
+ shape_v = [None, ] * self._dim
+ 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
# setup for rwork, iwork is useless.
if rwork is None:
# --- Local allocation ---
- self._rwork = [npw.zeros(shape_v)]
+ self._rwork = [npw.zeros(shape_v[d]) for d in xrange(self._dim)]
else:
assert isinstance(rwork, list), 'rwork must be a list.'
# --- External rwork ---
self._rwork = rwork
- assert len(self._rwork) == 1
- assert self._rwork[0].shape == shape_v
+ assert len(self._rwork) == self._dim
+ for i in xrange(self._dim):
+ assert self._rwork[i].shape == shape_v[i]
def _noca(self, dt):
"""
@@ -354,21 +344,21 @@ class NocaForces(Forces):
# Last part
# 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)
+ self._laplacian.fd_scheme.computeIndices_reduced(sl, buff.shape)
for j in i_t:
- self._rwork[...] = self._laplacian(vd[j], self._rwork)
+ buff[...] = self._laplacian(vd[j], buff)
res[i_n] += self._coeff * self.nu * normal\
- * npw.real_sum(coords[j] * self._rwork[sl])
- res[j] -= self._coeff * self.nu * normal * coords[i_n] * \
- npw.real_sum(self._rwork[sl])
- self._fd_scheme.computeIndices(sl)
- self._fd_scheme.compute(vd[i_n], i_n, self._rwork)
- res[i_n] += 2.0 * normal * self.nu * npw.real_sum(self._rwork[sl])
+ * npw.real_sum(coords[j] * buff)
+ res[j] -= self._coeff * self.nu * normal * x0 * \
+ npw.real_sum(buff)
+ self._fd_scheme.computeIndices_reduced(sl, buff.shape)
+ self._fd_scheme.compute(vd[i_n], i_n, buff)
+ res[i_n] += 2.0 * normal * self.nu * npw.real_sum(buff)
for j in i_t:
- self._fd_scheme.compute(vd[i_n], j, self._rwork)
- res[j] += normal * self.nu * npw.real_sum(self._rwork[sl])
- self._fd_scheme.compute(vd[j], i_n, self._rwork)
- res[j] += normal * self.nu * npw.real_sum(self._rwork[sl])
+ self._fd_scheme.compute(vd[i_n], j, buff)
+ res[j] += normal * self.nu * npw.real_sum(buff)
+ self._fd_scheme.compute(vd[j], i_n, buff)
+ res[j] += normal * self.nu * npw.real_sum(buff)
res *= dsurf
return res
diff --git a/HySoP/hysop/operator/drag_and_lift.py b/HySoP/hysop/operator/drag_and_lift.py
index 99ab378dca46f3859b8802f3f9655b1cbe6ff62a..2f97c4b9b8c2d6c1e93ed9aa54bb8c6cf454f835 100644
--- a/HySoP/hysop/operator/drag_and_lift.py
+++ b/HySoP/hysop/operator/drag_and_lift.py
@@ -37,7 +37,7 @@ class Forces(Computational):
vd = self.discreteFields[self.input[0]]
v_ind = vd.topology.mesh.iCompute
shape_v = vd[0][v_ind].shape
- return {'rwork': [shape_v], 'iwork': None}
+ return {'rwork': None, 'iwork': None}
def discretize(self):
super(Forces, self)._standard_discretize(self._min_ghosts)
@@ -159,10 +159,15 @@ class NocaForces(Forces):
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': [shape_v], 'iwork': None}
+
+ shape_v = [None, ] * self.domain.dimension
+ 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
+ # setup for rwork, iwork is useless.
+ return {'rwork': shape_v, 'iwork': None}
@opsetup
def setup(self, rwork=None, iwork=None):
diff --git a/HySoP/hysop/operator/tests/test_drag_and_lift.py b/HySoP/hysop/operator/tests/test_drag_and_lift.py
index 150c1019934fd5c90d9292d0517e1be2a4da35ea..c7645fa2e7570767d56f55f3b7709f7132726acd 100644
--- a/HySoP/hysop/operator/tests/test_drag_and_lift.py
+++ b/HySoP/hysop/operator/tests/test_drag_and_lift.py
@@ -75,6 +75,7 @@ def init_vort(discr, fileref):
vdref = veloref.discretize(topo)
return topo, wdref, vdref
+
def check_drag(penal, wref, vref, vorti, velo):
penal.discretize()
penal.setup()
@@ -109,7 +110,8 @@ def test_momentum_forces_3D():
penal.discretize()
penal.setup()
- dg = MomentumForces(velocity=velo, discretization=topo, subsets=[hsphere])
+ dg = MomentumForces(velocity=velo, discretization=topo,
+ obstacles=[hsphere])
dg.discretize()
dg.setup()
simu = Simulation(nbIter=3)