diff --git a/HySoP/hysop/domain/subsets/cylinder.py b/HySoP/hysop/domain/subsets/cylinder.py
index 1bcff7ef895e59eb561125b4ca60fdbd820cf795..a4dd309c9a152f8658a0fa450f7c17fa76985df5 100644
--- a/HySoP/hysop/domain/subsets/cylinder.py
+++ b/HySoP/hysop/domain/subsets/cylinder.py
@@ -54,15 +54,18 @@ class Cylinder(Subset):
             mesh = self._box.mesh[topo]
             shift = mesh.local_start
             dim = self._parent.dimension
-            tmp = np.where(np.logical_and(self.chi(*mesh.coords) <= 0,
-                                          mesh.coords[self.axis] ==
-                                          mesh.coords[self.axis]))
+            tmp = np.where(self.chi(*mesh.coords))
             self.ind[topo] = tuple([tmp[d] + shift[d] for d in xrange(dim)])
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
         return self.ind[topo]
 
+    def chi(self, *args):
+        return np.logical_and(self._is_inside(*args) <= 0,
+                              args[self.axis] ==
+                              args[self.axis])
+
     def __str__(self):
         s = self.__class__.__name__ + ' of radius ' + str(self.radius)
         s += ' and center position ' + str(self.origin)
@@ -90,16 +93,19 @@ class HemiCylinder(Cylinder):
             return x - cutpos
         self.LeftBox = LeftBox
 
+    def chi(self, *args):
+        return (np.logical_and(
+            np.logical_and(self._is_inside(*args) <= 0,
+                           self.LeftBox(args[self.cutdir]) <= 0),
+            args[self.axis] == args[self.axis]))
+
     def discretize(self, topo):
         if topo not in self.ind:
             self._box.discretize(topo)
             mesh = self._box.mesh[topo]
             shift = mesh.local_start
             dim = self._parent.dimension
-            tmp = np.where(np.logical_and(
-                np.logical_and(self.chi(*mesh.coords) <= 0,
-                               self.LeftBox(mesh.coords[self.cutdir]) <= 0),
-                mesh.coords[self.axis] == mesh.coords[self.axis]))
+            tmp = np.where(self.chi(*mesh.coords))
             self.ind[topo] = tuple([tmp[d] + shift[d] for d in xrange(dim)])
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
diff --git a/HySoP/hysop/domain/subsets/sphere.py b/HySoP/hysop/domain/subsets/sphere.py
index 5b5ea74dd1e91ad95dba06235e2c1c7b36f130d5..469b098546ced68fc202ea17c18b1467e473af81 100644
--- a/HySoP/hysop/domain/subsets/sphere.py
+++ b/HySoP/hysop/domain/subsets/sphere.py
@@ -44,7 +44,7 @@ class Sphere(Subset):
             mesh = self._box.mesh[topo]
             shift = mesh.local_start
             dim = self._parent.dimension
-            tmp = np.where(self.chi(*mesh.coords) <= 0)
+            tmp = np.where(self.chi(*mesh.coords))
             self.ind[topo] = tuple([tmp[d] + shift[d] for d in xrange(dim)])
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
@@ -79,6 +79,11 @@ class HemiSphere(Sphere):
         def LeftBox(x):
             return x - cutpos
         self.LeftBox = LeftBox
+        #self.chi = self._inside
+
+    def chi(self, *args):
+        return np.logical_and(self._is_inside(*args) <= 0,
+                              self.LeftBox(args[self.cutdir]) <= 0)
 
     def discretize(self, topo):
         assert isinstance(topo, Cartesian)
@@ -87,9 +92,9 @@ class HemiSphere(Sphere):
             mesh = self._box.mesh[topo]
             dim = self._parent.dimension
             shift = mesh.local_start
-            self.ind[topo] = np.where(
-                np.logical_and(self.chi(*mesh.coords) <= 0,
-                               self.LeftBox(mesh.coords[self.cutdir]) <= 0))
+            self.ind[topo] = np.where(self.chi(*mesh.coords))
+#                np.logical_and(self.chi(*mesh.coords) <= 0,
+#                               self.LeftBox(mesh.coords[self.cutdir]) <= 0))
             self.ind[topo] = tuple([self.ind[topo][d] + shift[d]
                                     for d in xrange(dim)])
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
diff --git a/HySoP/hysop/domain/subsets/subset.py b/HySoP/hysop/domain/subsets/subset.py
index f0102604158373ef05ed06706187e33569a29862..f22d1600ac2e8f3012b6e4cea794873b8b2a2797 100644
--- a/HySoP/hysop/domain/subsets/subset.py
+++ b/HySoP/hysop/domain/subsets/subset.py
@@ -44,7 +44,7 @@ class Subset(object):
     def __init__(self, parent, chi=None):
         """
         @param parent : the domain in which the subset is defined
-        @param chi : indicator function of the domain.
+        @param func : indicator function of the domain.
         """
         assert isinstance(parent, Domain)
         self._parent = parent
@@ -54,10 +54,13 @@ class Subset(object):
         ## Dictionnary (key = topo), on_proc[topo] = True
         ## if the subset has points on the current mpi process.
         self.on_proc = {}
-        ## indicator function of the subset
-        self.chi = chi
+        # indicator function of the subset
+        self._is_inside = chi
         self.is_porous = False
 
+    def chi(self, *args):
+        return self._is_inside(*args) <= 0
+
     def discretize(self, topo):
         """
         Create the list of indices for points inside the subset
@@ -67,7 +70,7 @@ class Subset(object):
         assert isinstance(topo, Cartesian)
         if topo not in self.ind:
             dim = self._parent.dimension
-            self.ind[topo] = np.where(self.chi(*topo.mesh.coords) <= 0)
+            self.ind[topo] = np.where(self.chi(*topo.mesh.coords))
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
@@ -91,13 +94,13 @@ class Subset(object):
         if len(slist) == 1:
             return sref.ind[topo]
         if sref.chi is not None:
-            c0 = sref.chi(*topo.mesh.coords) <= 0
+            c0 = sref.chi(*topo.mesh.coords)
         else:
             c0 = sref.mesh[topo].chi(*topo.mesh.coords)
 
         for s in slist[1:]:
             if s.chi is not None:
-                c1 = s.chi(*topo.mesh.coords) <= 0
+                c1 = s.chi(*topo.mesh.coords)
             else:
                 c1 = s.mesh[topo].chi(*topo.mesh.coords)
 
@@ -112,22 +115,32 @@ class Subset(object):
         @return : the union of a set of subsets for a
         given topo as a tuple of arrays, like self.ind[topo]
         """
+        return np.where(Subset.union_as_bool(slist, topo))
+
+    @staticmethod
+    def union_as_bool(slist, topo):
+        """
+        @param slist : a list of subsets
+        @param topo : a parmepy.mpi.topology.Cartesian
+        @return : the union of a set of subsets for a
+        given topo as an array of boolean
+        """
         sref = slist[0]
-        if len(slist) == 1:
-            return sref.ind[topo]
         if sref.chi is not None:
-            c0 = sref.chi(*topo.mesh.coords) <= 0
+            c0 = sref.chi(*topo.mesh.coords)
         else:
             c0 = sref.mesh[topo].chi(*topo.mesh.coords)
+        if len(slist) == 1:
+            return c0
 
         for s in slist[1:]:
             if s.chi is not None:
-                c1 = s.chi(*topo.mesh.coords) <= 0
+                c1 = s.chi(*topo.mesh.coords)
             else:
                 c1 = s.mesh[topo].chi(*topo.mesh.coords)
 
             c0 = np.logical_or(c0, c1)
-        return np.where(c0)
+        return c0
 
     @staticmethod
     def substraction(s1, s2, topo):
@@ -139,11 +152,25 @@ class Subset(object):
         given topo and return a tuple of arrays, like self.ind[topo]
         """
         if s1.chi is not None:
-            c1 = s1.chi(*topo.mesh.coords) <= 0
+            c1 = s1.chi(*topo.mesh.coords)
         else:
             c1 = s1.mesh[topo].chi(*topo.mesh.coords)
         if s2.chi is not None:
-            c2 = s2.chi(*topo.mesh.coords) > 0
+            c2 = np.logical_not(s2.chi(*topo.mesh.coords))
         else:
-            c2 = s2.mesh[topo].chi(*topo.mesh.coords) == False
+            c2 = np.logical_not(s2.mesh[topo].chi(*topo.mesh.coords))
+        return np.where(np.logical_and(c1, c2))
+
+    @staticmethod
+    def substract_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
+        other subsets for a
+        given topo and return a tuple of arrays, like self.ind[topo]
+        """
+        c1 = Subset.union_as_bool(s1, topo)
+        c2 = np.logical_not(Subset.union_as_bool(s2, topo))
         return np.where(np.logical_and(c1, c2))
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_2d.xmf b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_2d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..f2ac1aaee9384030068bfa2a638afc9869a22e08
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_2d.xmf
@@ -0,0 +1,25 @@
+<?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="2DCORECTMesh" NumberOfElements="96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDY">
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_2d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_2d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_2d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..78cea29e189ae798273add613c75b8430dde7d98
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_2d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_3d.xmf b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_3d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..5da7b0be0858505cc1aac470624ecaee9ba4d848
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_3d.xmf
@@ -0,0 +1,25 @@
+<?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="102  96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDYDZ">
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.5  -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.061599855952741041  0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="102  96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_3d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_3d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_3d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..a2be19b01d40c4b8aed359bc72a1bc50d5768530
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_3d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_2d.xmf b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_2d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..466751e4332b1257597a5bd969fd92ecf5c0c83f
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_2d.xmf
@@ -0,0 +1,25 @@
+<?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="2DCORECTMesh" NumberOfElements="96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDY">
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subs_2d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_2d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_2d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..3d748d574c75ca3e16a5b74784ce39b0d1d7ba82
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_2d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_3d.xmf b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_3d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..6070488b4adc923aa8186fa8d2b711b31f021e2e
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_3d.xmf
@@ -0,0 +1,25 @@
+<?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="102  96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDYDZ">
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.5  -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.061599855952741041  0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="102  96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subs_3d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_3d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_3d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..7471b641dc68c53f175e238cac3fe906ad8f3e86
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subs_3d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_2d.xmf b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_2d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..6a3d717135297b94a689f20f9d998bbae0f6186a
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_2d.xmf
@@ -0,0 +1,25 @@
+<?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="2DCORECTMesh" NumberOfElements="96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDY">
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subslist_2d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_2d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_2d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..a3f917d33c3d6894ba640a311218b26a255621a7
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_2d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_3d.xmf b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_3d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..ffff8875c81b06a9d890ded6fb2b5cb84307521e
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_3d.xmf
@@ -0,0 +1,25 @@
+<?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="102  96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDYDZ">
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.5  -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.061599855952741041  0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="102  96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subslist_3d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_3d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_3d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..0a8e08127c821f25d43b7e014d82473cc0d8aa7b
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p1/multi_obst_subslist_3d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_2d.xmf b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_2d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..f2ac1aaee9384030068bfa2a638afc9869a22e08
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_2d.xmf
@@ -0,0 +1,25 @@
+<?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="2DCORECTMesh" NumberOfElements="96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDY">
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_2d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_2d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_2d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..62f1897e48c5d545321af4295bef85b8ff03eebd
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_2d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_3d.xmf b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_3d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..5da7b0be0858505cc1aac470624ecaee9ba4d848
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_3d.xmf
@@ -0,0 +1,25 @@
+<?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="102  96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDYDZ">
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.5  -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.061599855952741041  0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="102  96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_3d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_3d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_3d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..70bc9415d5088a2598588fb9bff5f2bff5cc161c
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_3d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_2d.xmf b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_2d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..466751e4332b1257597a5bd969fd92ecf5c0c83f
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_2d.xmf
@@ -0,0 +1,25 @@
+<?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="2DCORECTMesh" NumberOfElements="96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDY">
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subs_2d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_2d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_2d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..9f62457b64a5663ffbfe09a82b6cdf5d8f2dbe0c
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_2d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_3d.xmf b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_3d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..6070488b4adc923aa8186fa8d2b711b31f021e2e
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_3d.xmf
@@ -0,0 +1,25 @@
+<?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="102  96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDYDZ">
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.5  -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.061599855952741041  0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="102  96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subs_3d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_3d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_3d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..e86ec3a2fde4456cb9fdcc3b1d974172105868f3
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subs_3d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_2d.xmf b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_2d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..6a3d717135297b94a689f20f9d998bbae0f6186a
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_2d.xmf
@@ -0,0 +1,25 @@
+<?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="2DCORECTMesh" NumberOfElements="96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDY">
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="2 " NumberType="Float" Precision="8" Format="XML">
+     0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subslist_2d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_2d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_2d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..0a2b8be0d1ee09de08be01925d900a26ff908dcf
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_2d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_3d.xmf b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_3d.xmf
new file mode 100644
index 0000000000000000000000000000000000000000..ffff8875c81b06a9d890ded6fb2b5cb84307521e
--- /dev/null
+++ b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_3d.xmf
@@ -0,0 +1,25 @@
+<?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="102  96  128 "/>
+    <Geometry GeometryType="ORIGIN_DXDYDZ">
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.5  -0.29999999999999999  0.10000000000000001
+     </DataItem>
+     <DataItem Dimensions="3 " NumberType="Float" Precision="8" Format="XML">
+     0.061599855952741041  0.065449846949787352  0.049087385212340517
+     </DataItem>
+    </Geometry>
+    <Attribute Name="vref_0_X" AttributeType="Scalar" Center="Node">
+     <DataItem Dimensions="102  96  128 " NumberType="Float" Precision="8" Format="HDF" Compression="Raw">
+      multi_obst_subslist_3d_00000.h5:/vref_0_X
+     </DataItem>
+    </Attribute>
+   </Grid>
+  </Grid>
+ </Domain>
+</Xdmf>
diff --git a/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_3d_00000.h5 b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_3d_00000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..7dc3bf1c9148324d8b5f31aaca154dcc00020c82
Binary files /dev/null and b/HySoP/hysop/domain/tests/ref_files/p8/multi_obst_subslist_3d_00000.h5 differ
diff --git a/HySoP/hysop/domain/tests/test_subset.py b/HySoP/hysop/domain/tests/test_subset.py
index 8cc16476cc151e0ffbd7f9efe58b693d6c2f7d75..555916e67c961c5e3dfcbf5e6b7d228bb67c9d10 100644
--- a/HySoP/hysop/domain/tests/test_subset.py
+++ b/HySoP/hysop/domain/tests/test_subset.py
@@ -1,4 +1,6 @@
 from parmepy.domain.subsets.sphere import Sphere, HemiSphere
+from parmepy.domain.subsets.subset import Subset
+from parmepy.domain.subsets.boxes import SubBox
 from parmepy.domain.subsets.cylinder import Cylinder, HemiCylinder
 from parmepy.operator.hdf_io import HDF_Reader
 from parmepy.tools.parameters import Discretization, IO_params
@@ -36,7 +38,7 @@ def check_subset(discr, fileref, class_name):
     vref = Field(domain=topo.domain, name='vref')
     vd = vref.discretize(topo)
     iop = IO_params(working_dir + fileref)
-    reader = HDF_Reader(variables={vref: topo}, io_params=iop)
+    reader = HDF_Reader(variables={vref: topo}, io_params=iop, restart=0)
     reader.discretize()
     reader.setup()
     reader.apply()
@@ -53,31 +55,141 @@ def check_subset(discr, fileref, class_name):
     scal = Field(domain=topo.domain, isVector=False)
     sd = scal.discretize(topo)
     sd[0][subs.ind[topo]] = 2.
-    return np.allclose(sd[0], vd[0])
+    iCompute = topo.mesh.iCompute
+    return np.allclose(sd[0][iCompute], vd[0][iCompute])
 
 
 def test_sphere_3d():
-    assert check_subset(discr3D, 'sphere3d_00000.h5', Sphere)
+    assert check_subset(discr3D, 'sphere3d', Sphere)
 
 
 def test_hemisphere_3d():
-    assert check_subset(discr3D, 'hemisphere3d_00000.h5', HemiSphere)
+    assert check_subset(discr3D, 'hemisphere3d', HemiSphere)
 
 
 def test_cylinder_3d():
-    assert check_subset(discr3D, 'cylinder3d_00000.h5', Cylinder)
+    assert check_subset(discr3D, 'cylinder3d', Cylinder)
 
 
 def test_hemicylinder_3d():
-    assert check_subset(discr3D, 'hemicylinder3d_00000.h5', HemiCylinder)
+    assert check_subset(discr3D, 'hemicylinder3d', HemiCylinder)
 
 
 def test_sphere_2d():
-    assert check_subset(discr2D, 'sphere2d_00000.h5', Sphere)
+    assert check_subset(discr2D, 'sphere2d', Sphere)
 
 
 def test_hemisphere_2d():
-    assert check_subset(discr2D, 'hemisphere2d_00000.h5', HemiSphere)
+    assert check_subset(discr2D, 'hemisphere2d', HemiSphere)
+
+
+def init_multi(discr, fileref):
+    Cartesian.reset_counter()
+    dim = len(discr.resolution)
+    dom = Box(dimension=dim, length=ldom[:dim],
+              origin=xdom[:dim])
+    topo = dom.create_topology(discr)
+    vref = Field(domain=topo.domain, name='vref', isVector=False)
+    vd = vref.discretize(topo)
+    iop = IO_params(working_dir + fileref)
+    reader = HDF_Reader(variables={vref: topo}, io_params=iop, restart=0)
+    reader.discretize()
+    reader.setup()
+    reader.apply()
+    reader.finalize()
+    scal = Field(domain=topo.domain, isVector=False)
+    sd = scal.discretize(topo)
+    return topo, sd, vd
+
+
+def init_obst_list(discr, fileref):
+    topo, sd, vd = init_multi(discr, fileref)
+    rd = ldom[0] * 0.1
+    dim = len(discr.resolution)
+    s1 = Sphere(origin=xpos[:dim], radius=rd * 1.4, parent=topo.domain)
+    s2 = Cylinder(origin=xpos[:dim], radius=rd, parent=topo.domain)
+    newpos = xpos.copy()
+    newpos[1] += 3 * rd
+    s3 = Sphere(origin=newpos[:dim], radius=1.3 * rd, parent=topo.domain)
+    obs_list = [s1, s2, s3]
+    for obs in obs_list:
+        obs.discretize(topo)
+    sd = Field(domain=topo.domain, isVector=False).discretize(topo)
+    indices = Subset.union(obs_list, topo)
+    sd[0][indices] = 2.
+    iCompute = topo.mesh.iCompute
+    return np.allclose(sd[0][iCompute], vd[0][iCompute])
+
+
+def test_union_3D():
+    assert init_obst_list(discr3D, 'multi_obst_3d')
+
+
+def test_union_2D():
+    assert init_obst_list(discr2D, 'multi_obst_2d')
+
+
+def init_substract(discr, fileref):
+    topo, sd, vd = init_multi(discr, fileref)
+    dim = len(discr.resolution)
+    rd = ldom[0] * 0.1
+    s1 = Sphere(origin=xpos[:dim], radius=rd, parent=topo.domain)
+    pos = topo.domain.origin + 0.05
+    ll = topo.domain.length - 0.2
+    box = SubBox(origin=pos, length=ll, parent=topo.domain)
+    s1.discretize(topo)
+    box.discretize(topo)
+    sd = Field(domain=topo.domain, isVector=False).discretize(topo)
+    indices = Subset.substraction(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_substraction_2D():
+    assert init_substract(discr2D, 'multi_obst_subs_2d')
+
+
+def init_substract_lists(discr, fileref):
+    topo, sd, vd = init_multi(discr, fileref)
+    dim = len(discr.resolution)
+    rd = ldom[0] * 0.1
+    obs = []
+    obs.append(Sphere(origin=xpos[:dim], radius=rd, parent=topo.domain))
+    pos = topo.domain.origin + 0.05
+    ll = topo.domain.length - 0.2
+    box = []
+    box.append(SubBox(origin=pos, length=2 * ll / 3., parent=topo.domain))
+    pos2 = topo.domain.origin + 0.05
+    pos2[1] += topo.domain.length[1] / 5
+    box.append(SubBox(origin=pos2, length=4 * ll / 5., parent=topo.domain))
+    obs.append(Cylinder(origin=xpos[:dim], radius=rd * 0.5,
+                        parent=topo.domain))
+    newpos = xpos.copy()
+    newpos[1] += 3 * rd
+    obs.append(Sphere(origin=newpos[:dim], radius=1.3 * rd,
+                      parent=topo.domain))
+    for ob in obs:
+        ob.discretize(topo)
+    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)
+    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_substraction_list_2D():
+    assert init_substract_lists(discr2D, 'multi_obst_subslist_2d')
+
 
 # This may be useful to run mpi tests
 if __name__ == "__main__":
@@ -87,3 +199,9 @@ if __name__ == "__main__":
     test_hemicylinder_3d()
     test_sphere_2d()
     test_hemisphere_2d()
+    test_union_2D()
+    test_union_3D()
+    test_substraction_3D()
+    test_substraction_2D()
+    test_substraction_list_3D()
+    test_substraction_list_2D()
diff --git a/HySoP/hysop/f2py/scales2py.f90 b/HySoP/hysop/f2py/scales2py.f90
index 8171a02eb5bf6cc9b96d0ded3dfaa88e92bb83cd..78845354729494b91c0f97220010d54e4ab4dd06 100755
--- a/HySoP/hysop/f2py/scales2py.f90
+++ b/HySoP/hysop/f2py/scales2py.f90
@@ -35,7 +35,7 @@ contains
     !f2py intent(hide) dim
     !f2py character(*) optional, intent(in) :: order = 'p_O2'
     !f2py, depends(dim), intent(in) :: topodims
-    integer :: error,groupsize !rank, nbprocs
+    integer :: error  !,groupsize !rank, nbprocs
 
     if(dim /= 3) then
        stop 'Scales advection solver initialisation failed : not yet implemented for 2d problem.'
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index a22a48a758e11d5cb864976e6ab3ea2e04a629fc..ad272ceec21dce3a1270ee622099e3f4c7aef95c 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -20,21 +20,18 @@ class DragAndLift(DiscreteOperator):
     The present class implements formula (52) of Plouhmans2002.
     Integral inside the obstacle is not taken into account.
     """
-    def __init__(self, velocity, vorticity, nu, normalization,
-                 volumeOfControl, obstacles=None, **kwds):
+    def __init__(self, velocity, vorticity, nu, volumeOfControl,
+                 normalization=1., obstacles=None, **kwds):
         """
         @param velocity field
         @param vorticity field
         @param nu viscosity
-        @param the topology on which forces will be computed
         @param a volume of control
-        (parmepy.domain.obstacle.controlBox.ControlBox object)
-        @param obstacles a list of parmepy.domain.obstacles inside
-        the box
-        @param io_params : parameters (dict) to set file output.
-        If  None, no output. Set io_params = {} if you want output,
-        with default parameters values. Default file name = 'drag_and_lift'
-        See parmepy.tools.io_utils.Writer for details
+        @param normalization : a normalization coefficient applied to the force
+        (default = 1.)
+        (parmepy.domain.subset.boxes.SubBox)
+        @param obstacles a list of parmepy.domain.subsets inside the volume of
+        control
         """
         assert 'variables' not in kwds, 'variables parameter is useless.'
         super(DragAndLift, self).__init__(variables=[velocity, vorticity],
@@ -77,19 +74,43 @@ class DragAndLift(DiscreteOperator):
         # a scalar field in a given direction.
         # Default = FD_C_2. Todo : set this as an input method value.
         self._fd_scheme = FD_C_2(self._topology.mesh.space_step)
+
+        # Get the list of indices of points used to compute the
+        # integrals over the volume.
+        self._indices = self._init_indices(obstacles)
+
+        # prepare ghost points synchro for velocity and vorticity used
+        # in fd schemes
+        nbc = self.velocity.nbComponents + self.vorticity.nbComponents
+        self._synchronize = UpdateGhosts(self._topology, nbc)
+
+        # Which formula must be used to compute the forces ...
+        self._formula = self._switch_method()
+
+        self._on_proc = self._voc.on_proc[self._topology]
+
+    def _init_indices(self, obstacles):
+        """
+        Compute a list of indices corresponding to points inside
+        the volume of control minus those inside the obstacles
+        """
+        from parmepy.domain.subsets.subset import Subset
         self._voc.discretize(self._topology)
         if obstacles is None:
-            self._obsinds = []
+            return self._voc.ind[self._topology]
         else:
             for obs in obstacles:
                 obs.discretize(self._topology)
-            self._obsinds = [ind for ind in obs.ind[self._topology]
-                             for obs in obstacles]
+            return Subset.substract_list_of_sets([self._voc], obstacles,
+                                                 self._topology)
 
-        # prepare ghost points synchro for velocity and vorticity used
-        # in fd schemes
-        nbc = self.velocity.nbComponents + self.vorticity.nbComponents
-        self._synchronize = UpdateGhosts(self._topology, nbc)
+    def _switch_method(self):
+        """
+        select the proper method to compute drag and lift
+        """
+        # Heloise formula
+        # assert self._dim == 3, 'Not defined for dim < 3'
+        return self._momentum
 
     def _set_work_arrays(self, rwork=None, iwork=None):
 
@@ -123,7 +144,7 @@ class DragAndLift(DiscreteOperator):
 
     @debug
     @profile
-    def apply(self, simulation=None):
+    def _noca(self, simulation=None):
         """
         Perform integrals on volume and surfaces of the control box
         @param parmepy.problem.simulation : object describing
@@ -169,6 +190,32 @@ class DragAndLift(DiscreteOperator):
 
         return self.force
 
+    def _momentum(self, dt):
+        # Synchro of ghost points is required for fd schemes
+        self._synchronize(self.velocity.data + self.vorticity.data)
+
+        # -- Integration over the volume of control --
+        if not self._on_proc:
+            self._buffer[...] = 0.0
+            self.force[...] = 0.0
+        else:
+            self._buffer = self._voc.integrate_dfield_allc(self.velocity)
+            self.force[...] = 1. / dt * (self._buffer - self._previous)
+
+        # Update previous for next time step ...
+        self._previous[...] = self._buffer[...]
+
+    def apply(self, simulation=None):
+        """
+        Perform integrals on volume and surfaces of the control box
+        @param parmepy.problem.simulation : object describing
+        simulation parameters
+        """
+        assert simulation is not None,\
+            "Simulation parameter is required for DragAndLift apply."
+        dt = simulation.timeStep
+        self._formula(dt)
+
     def _integrateOnSurface(self, surf, res):
 
         res[...] = 0.0
diff --git a/HySoP/hysop/operator/drag_and_lift.py b/HySoP/hysop/operator/drag_and_lift.py
index c92be7f60822b27850157a58e69b98af21e88360..4d72c35c34811d015381b732900c2c7963555704 100644
--- a/HySoP/hysop/operator/drag_and_lift.py
+++ b/HySoP/hysop/operator/drag_and_lift.py
@@ -14,8 +14,8 @@ class DragAndLift(Computational):
     The present class implements formula (52) of Plouhmans2002.
     Integral inside the obstacle is not taken into account.
     """
-    def __init__(self, velocity, vorticity, nu, normalization,
-                 volumeOfControl, obstacles=None, **kwds):
+    def __init__(self, velocity, vorticity, volumeOfControl,
+                 nu=0., normalization=1., obstacles=None, **kwds):
         """
         @param velocity field
         @param vorticity field
diff --git a/HySoP/hysop/operator/tests/ref_files/p1/penal2d_multi_00000.h5 b/HySoP/hysop/operator/tests/ref_files/p1/penal2d_multi_00000.h5
index c4db3720d9290e32279006f0fbf42ee0e042bf8f..f1f8b63899a516bf5136e4c3d9d36d7aa7af0ac5 100644
Binary files a/HySoP/hysop/operator/tests/ref_files/p1/penal2d_multi_00000.h5 and b/HySoP/hysop/operator/tests/ref_files/p1/penal2d_multi_00000.h5 differ
diff --git a/HySoP/hysop/operator/tests/ref_files/p1/penal3d_multi_00000.h5 b/HySoP/hysop/operator/tests/ref_files/p1/penal3d_multi_00000.h5
index bd3c49c0867bd83014a2dc9c8404fd768338be02..3701ce5ba5e2bb7738660b8ab6fc975269afb307 100644
Binary files a/HySoP/hysop/operator/tests/ref_files/p1/penal3d_multi_00000.h5 and b/HySoP/hysop/operator/tests/ref_files/p1/penal3d_multi_00000.h5 differ
diff --git a/HySoP/hysop/operator/tests/ref_files/p1/penal3d_porous_cyl_00000.h5 b/HySoP/hysop/operator/tests/ref_files/p1/penal3d_porous_cyl_00000.h5
index 7826d569c34a57e5f59e8b9cb4470f251ea54d73..a60f349e5ed5925bae9466e9763bd39ef32f9b3f 100644
Binary files a/HySoP/hysop/operator/tests/ref_files/p1/penal3d_porous_cyl_00000.h5 and b/HySoP/hysop/operator/tests/ref_files/p1/penal3d_porous_cyl_00000.h5 differ
diff --git a/HySoP/hysop/operator/tests/ref_files/p8/penal2d_multi_00000.h5 b/HySoP/hysop/operator/tests/ref_files/p8/penal2d_multi_00000.h5
index 1b0f74d6487592c2515a3595cb61ada163414833..85c94635f570c6da50d12b41df46b06ccb5fcfed 100644
Binary files a/HySoP/hysop/operator/tests/ref_files/p8/penal2d_multi_00000.h5 and b/HySoP/hysop/operator/tests/ref_files/p8/penal2d_multi_00000.h5 differ
diff --git a/HySoP/hysop/operator/tests/ref_files/p8/penal3d_multi_00000.h5 b/HySoP/hysop/operator/tests/ref_files/p8/penal3d_multi_00000.h5
index 251473bcd78973d7bf0e1c4635429d1130a896ff..5347f06b0464fd59c1b53a0e3e01eb2fead1d046 100644
Binary files a/HySoP/hysop/operator/tests/ref_files/p8/penal3d_multi_00000.h5 and b/HySoP/hysop/operator/tests/ref_files/p8/penal3d_multi_00000.h5 differ
diff --git a/HySoP/hysop/operator/tests/ref_files/p8/penal3d_porous_cyl_00000.h5 b/HySoP/hysop/operator/tests/ref_files/p8/penal3d_porous_cyl_00000.h5
index d78632277ffb1cccdf1d5e8f8d8eb5a4efd7e71f..af3761b545ab8b901170664d6dc48e7e1c6a67a2 100644
Binary files a/HySoP/hysop/operator/tests/ref_files/p8/penal3d_porous_cyl_00000.h5 and b/HySoP/hysop/operator/tests/ref_files/p8/penal3d_porous_cyl_00000.h5 differ
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 148dae2686d67df2f265964223174ca3c6724c89..f61a0d3c413631602ab7a1b6f3a5e1b562bd08c4 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -115,8 +115,8 @@ if("@WITH_GPU@" is "ON"):
 
 config = Configuration(name=name,
       version='@PYPACKAGE_VERSION@',
-      description='Particular Methods implementation in Python',
-      author='Jean-Matthieu Etancelin, Franck Pérignon, Christophe Picard',
+      description='Particular Methods implementation for heterogenous platforms.',
+      author='G.H Cottet, J.M Etancelin, C.Mimeau, F.Pérignon, C. Picard',
       author_email='parmes-devel@lists.forge.imag.fr',
       url='https://forge.imag.fr/projects/parmes/',
       license='GNU public license',