From 0249d5bf3853a0aa002a75a97b8408e899b29af4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 20 May 2016 17:01:52 +0200
Subject: [PATCH] Fix some tests (failed since previous commit) + pep8 things

---
 examples/NSDebug_faux2D.py                    |  2 +-
 examples/testCurl.py                          |  4 +-
 examples/testVisu.py                          |  2 +-
 hysop/domain/mesh.py                          | 45 +++++++-----
 hysop/domain/submesh.py                       | 16 ++---
 hysop/domain/subsets.py                       |  2 +-
 hysop/domain/tests/test_mesh.py               |  6 +-
 hysop/domain/tests/test_porous.py             | 10 +--
 hysop/domain/tests/test_regular_subset.py     |  2 +-
 hysop/domain/tests/test_submesh.py            |  2 +-
 hysop/domain/tests/test_subset.py             |  8 +--
 hysop/fields/discrete.py                      |  4 +-
 .../tests/test_gpu_multiresolution_filter.py  | 32 ++++-----
 hysop/mpi/bridge_inter.py                     |  4 +-
 hysop/mpi/topology.py                         |  4 +-
 hysop/numerics/differential_operations.py     |  6 +-
 hysop/numerics/finite_differences.py          |  4 +-
 hysop/numerics/tests/test_diffOp.py           |  4 +-
 hysop/numerics/tests/test_update_ghosts.py    |  2 +-
 hysop/numerics/utils.py                       |  2 +-
 hysop/operator/discrete/baroclinic.py         | 68 +++++++++----------
 .../operator/discrete/baroclinic_from_rhs.py  |  4 +-
 hysop/operator/discrete/density.py            |  8 +--
 hysop/operator/discrete/energy_enstrophy.py   | 12 ++--
 hysop/operator/discrete/monitoringPoints.py   | 10 +--
 hysop/operator/discrete/multiphase_gradp.py   | 34 +++++-----
 .../discrete/multiresolution_filter.py        | 12 ++--
 hysop/operator/discrete/profiles.py           | 10 +--
 hysop/operator/discrete/residual.py           | 10 +--
 .../operator/discrete/velocity_correction.py  |  4 +-
 hysop/operator/energy_enstrophy.py            |  4 +-
 hysop/operator/hdf_io.py                      |  2 +-
 hysop/operator/monitoringPoints.py            |  4 +-
 hysop/operator/profiles.py                    |  4 +-
 hysop/operator/residual.py                    |  2 +-
 hysop/operator/tests/test_differential.py     |  2 +-
 hysop/operator/tests/test_hdf5_io.py          | 16 ++---
 hysop/operator/tests/test_multiphase_gradp.py |  8 +--
 .../tests/test_multiresolution_filter.py      | 26 +++----
 hysop/operator/tests/test_penalization.py     |  4 +-
 hysop/operator/tests/test_redistribute.py     | 18 ++---
 41 files changed, 218 insertions(+), 205 deletions(-)

diff --git a/examples/NSDebug_faux2D.py b/examples/NSDebug_faux2D.py
index a6540feed..16902a658 100755
--- a/examples/NSDebug_faux2D.py
+++ b/examples/NSDebug_faux2D.py
@@ -373,7 +373,7 @@ cb2 = op['correction'].cb
 def pseudo_norm(topo, var):
     sloc = npw.zeros(3)
     gloc = npw.zeros(3)
-    sl = topo.mesh.iCompute
+    sl = topo.mesh.compute_index
     for i in xrange(3):
         sloc[i] = np.sum((var[i][sl]))
         # sfft[i] = la.norm((vfft[i][slfft]))
diff --git a/examples/testCurl.py b/examples/testCurl.py
index 7b0b93cc9..0a08866ab 100644
--- a/examples/testCurl.py
+++ b/examples/testCurl.py
@@ -98,8 +98,8 @@ curlFD.apply(simulation)
 printer.apply(simulation)
 printer2.apply(simulation)
 
-ind = topo3.mesh.iCompute
-indfft = topofft.mesh.iCompute
+ind = topo3.mesh.compute_index
+indfft = topofft.mesh.compute_index
 
 
 for i in xrange(3):
diff --git a/examples/testVisu.py b/examples/testVisu.py
index e2c7504a4..a7b10a5e1 100755
--- a/examples/testVisu.py
+++ b/examples/testVisu.py
@@ -125,7 +125,7 @@ for pref in filelist[1:]:
     simu.advance()
 ##     #print scalRef.norm(topo)
 ##     #print scal3D.norm(topo)
-##     #assert np.allclose(sc3D.data[0][topo.mesh.iCompute], scRef.data[0][topo.mesh.iCompute])
+##     #assert np.allclose(sc3D.data[0][topo.mesh.compute_index], scRef.data[0][topo.mesh.compute_index])
 
 
 velo = Field(domain=dom, name='Velocity', is_vector=True)
diff --git a/hysop/domain/mesh.py b/hysop/domain/mesh.py
index 4b0ec7576..01ea25d86 100644
--- a/hysop/domain/mesh.py
+++ b/hysop/domain/mesh.py
@@ -111,8 +111,8 @@ class Mesh(object):
         # the number of ghost points on each subdomain.
         self.discretization = discretization
         # Mesh step size
-        self.space_step = self.domain.length / (self.discretization.resolution
-                                                - 1)
+        self.space_step = self.domain.length / (
+            self.discretization.resolution - 1)
         # Coordinates of the origin of the mesh
         self.global_origin = self.domain.origin
         # Length of this mesh
@@ -133,6 +133,10 @@ class Mesh(object):
         start = npw.asdimarray(global_start)
         # last computational point, global index
         ghosts = self.discretization.ghosts
+        # resolution of the computational grid (i.e. grid without ghosts)
+        self.compute_resolution = tuple(npw.const_dimarray(resolution) -
+                                        2 * ghosts)
+
         # shift between local and global index,
         # such that: iglobal + shift = ilocal.
         # Usefull for start() and end() methods.
@@ -149,16 +153,20 @@ class Mesh(object):
         start = self.discretization.ghosts
         # last computational point, local index
         stop = self.resolution - start
-        # Local position of the computational points
-        self.iCompute = [slice(start[d], stop[d])
-                         for d in xrange(self._dim)]
+        # Local indices of the computational points
+        self.compute_index = [slice(start[d], stop[d])
+                              for d in xrange(self._dim)]
+        # Local indices of all points
+        self.local_index = [slice(0, self.resolution[d])
+                            for d in xrange(self._dim)]
 
         # coordinates of the points of the grid
         self.coords = np.ix_(*(np.linspace(self.origin[d],
                                            self.end[d], self.resolution[d])
                                for d in xrange(self._dim)))
         # coordinates of the 'computational' points, i.e. excluding ghosts.
-        self.compute_coords = self.reduce_coords(self.coords, self.iCompute)
+        self.compute_coords = self.reduce_coords(self.coords,
+                                                 self.compute_index)
 
         gend = self.discretization.resolution - 1
         # Is this mesh on the last process in some direction in the
@@ -166,7 +174,7 @@ class Mesh(object):
         is_last = np.asarray([gend[d] < self.position[d].stop
                               for d in xrange(self._dim)])
         # indices of points that must be used to compute integrals on this mesh
-        self.ind4integ = self.compute_integ_point(is_last, self.iCompute)
+        self.ind4integ = self.compute_integ_point(is_last, self.compute_index)
         self.coords4int = self.reduce_coords(self.coords, self.ind4integ)
         # True if the current mesh is locally defined on the current mpi proc.
         self.on_proc = True
@@ -195,8 +203,7 @@ class Mesh(object):
         if n_dir is not None:
             stops[n_dir] = npw.asdimarray([ic[d].start + 1 for d in n_dir])
 
-        return [slice(ic[d].start,
-                stops[d]) for d in xrange(dim)]
+        return [slice(ic[d].start, stops[d]) for d in xrange(dim)]
 
     @staticmethod
     def reduce_coords(coords, reduced_index):
@@ -227,20 +234,24 @@ class Mesh(object):
 
     def local_indices(self, point):
         """
-        returns indices of the point of coordinates (close to) tab = x, y, z
+        returns indices of the point of coordinates (close to) point = x, y, z
         If (x, y, z) is not a grid point, it returns the closest grid point.
         """
         point = npw.asrealarray(point)
         if self.is_inside(point):
-            return npw.asdimarray(np.rint((point - self.origin)
-                                          / self.space_step))
+            return npw.asdimarray(np.rint((point - self.origin) /
+                                          self.space_step))
         else:
             return False
 
     def global_indices(self, point):
+        """returns global indices of the point of coordinates
+        (close to) point = x, y, z
+        If (x, y, z) is not a grid point, it returns the closest grid point.
+        """
         point = npw.asrealarray(point)
-        return npw.asdimarray(np.rint((point - self.domain.origin)
-                                      / self.space_step))
+        return npw.asdimarray(np.rint((point - self.domain.origin) /
+                                      self.space_step))
 
     def is_inside(self, point):
         """True if the point belongs to volume or surface
@@ -261,7 +272,7 @@ class Mesh(object):
         s += ' - resolution : ' + str(self.resolution) + '\n'
         s += ' - position in global mesh: ' + str(self.position) + '\n'
         s += ' - global discretization : ' + str(self.discretization) + '\n'
-        s += ' - computation points :' + str(self.iCompute) + '\n'
+        s += ' - computation points :' + str(self.compute_index) + '\n'
         return s
 
     def convert2local(self, sl):
@@ -334,5 +345,7 @@ class Mesh(object):
         return self.discretization.resolution - 1
 
     def local_shift(self, indices):
-        shift = [self.iCompute[d].start for d in xrange(self._dim)]
+        """Shift input indices with ghost layer size
+        """
+        shift = [self.compute_index[d].start for d in xrange(self._dim)]
         return tuple([indices[d] + shift[d] for d in xrange(self._dim)])
diff --git a/hysop/domain/submesh.py b/hysop/domain/submesh.py
index 46b0735fe..412f5fb63 100644
--- a/hysop/domain/submesh.py
+++ b/hysop/domain/submesh.py
@@ -88,16 +88,16 @@ class SubMesh(object):
             self.position_in_parent = [s for s in sl]
             # Indices of the points of the submesh, relative to
             # the LOCAL array
-            self.iCompute = self.mesh.convert2local(self.position_in_parent)
+            self.compute_index = self.mesh.convert2local(self.position_in_parent)
 
             # Resolution of the local submesh
-            self.resolution = [self.iCompute[d].stop - self.iCompute[d].start
+            self.resolution = [self.compute_index[d].stop - self.compute_index[d].start
                                for d in xrange(self._dim)]
 
-            # Same as self.iCompute but recomputed to be used
+            # Same as self.compute_index but recomputed to be used
             # for integration on the submesh
             self.ind4integ = self.mesh.compute_integ_point(is_last,
-                                                           self.iCompute,
+                                                           self.compute_index,
                                                            self._n_dir)
             start = [self.position_in_parent[d].start - self.substart[d]
                      for d in xrange(self._dim)]
@@ -108,16 +108,16 @@ class SubMesh(object):
         else:
             self.position_in_parent = None
             self.position = None
-            self.iCompute = [slice(0, 0), ] * self._dim
-            self.ind4integ = self.iCompute
+            self.compute_index = [slice(0, 0), ] * self._dim
+            self.ind4integ = self.compute_index
             self.resolution = [0, ] * self._dim
 
         # Shift between local submesh and local parent mesh.
-        self.local_start = [self.iCompute[d].start for d in xrange(self._dim)]
+        self.local_start = [self.compute_index[d].start for d in xrange(self._dim)]
 
         # Coordinates of the points of the local submesh
         self.coords = self.mesh.reduce_coords(self.mesh.coords,
-                                              self.iCompute)
+                                              self.compute_index)
         self.coords4int = self.mesh.reduce_coords(self.mesh.coords,
                                                   self.ind4integ)
 
diff --git a/hysop/domain/subsets.py b/hysop/domain/subsets.py
index 02a36beb2..5dbdd9b24 100644
--- a/hysop/domain/subsets.py
+++ b/hysop/domain/subsets.py
@@ -744,7 +744,7 @@ class SubBox(Subset):
         self.mesh[topo] = SubMesh(topo.mesh, gstart, gend)
         self.real_length[topo] = self.mesh[topo].global_length
         self.real_orig[topo] = self.mesh[topo].global_origin
-        self.ind[topo] = [self.mesh[topo].iCompute]
+        self.ind[topo] = [self.mesh[topo].compute_index]
         self._reduce_topology(topo)
         return self.ind[topo]
 
diff --git a/hysop/domain/tests/test_mesh.py b/hysop/domain/tests/test_mesh.py
index 0ae2b4868..261abb444 100755
--- a/hysop/domain/tests/test_mesh.py
+++ b/hysop/domain/tests/test_mesh.py
@@ -35,7 +35,7 @@ def create_mesh(discr):
     # Test local grid
     assert (mesh.resolution == resolution).all()
     assert (mesh.origin == dom.origin + (shift - gh) * mesh.space_step).all()
-    assert mesh.iCompute == [slice(gh[d], resolution[d] - gh[d])
+    assert mesh.compute_index == [slice(gh[d], resolution[d] - gh[d])
                              for d in xrange(dim)]
     assert (mesh.global_indices(dom.origin) == [0, ] * dim).all()
     assert (mesh.local_indices(mesh.origin) == [0, ] * dim).all()
@@ -97,9 +97,9 @@ def test_convert_local():
     elif main_size == 8:
         newstart = start - topo.rank * Nz / main_size
         newend = end - topo.rank * Nz / main_size
-        slref = slice(max(mesh.iCompute[-1].start, newstart + gh[-1]),
+        slref = slice(max(mesh.compute_index[-1].start, newstart + gh[-1]),
                       min(max(newend + gh[-1], gh[-1]),
-                          mesh.iCompute[-1].stop))
+                          mesh.compute_index[-1].stop))
         sl = [slice(start + gh[d], end + gh[d]) for d in xrange(dim)]
         sl[-1] = slref
         if res is not None:
diff --git a/hysop/domain/tests/test_porous.py b/hysop/domain/tests/test_porous.py
index 3ae7c79c9..39dd869dd 100755
--- a/hysop/domain/tests/test_porous.py
+++ b/hysop/domain/tests/test_porous.py
@@ -54,7 +54,7 @@ def assert_porous(dref, sourcename, filename):
     for i in xrange(li):
         sd.data[0][ind[i]] = 10 * i + 1
     sdref = sref.discretize(topo)
-    ic = topo.mesh.iCompute
+    ic = topo.mesh.compute_index
     assert np.allclose(sd.data[0][ic], sdref.data[0][ic])
 
 
@@ -95,7 +95,7 @@ def assert_bipole(dref, sourcename, filename):
     for i in xrange(li):
         sd.data[0][ind[i]] = 10 * i + 1
     sdref = sref.discretize(topo)
-    ic = topo.mesh.iCompute
+    ic = topo.mesh.compute_index
     assert np.allclose(sd.data[0][ic], sdref.data[0][ic])
 
 
@@ -136,7 +136,7 @@ def assert_quadripole(dref, sourcename, filename):
     for i in xrange(li):
         sd.data[0][ind[i]] = 10 * i + 1
     sdref = sref.discretize(topo)
-    ic = topo.mesh.iCompute
+    ic = topo.mesh.compute_index
     assert np.allclose(sd.data[0][ic], sdref.data[0][ic])
 
 
@@ -161,7 +161,7 @@ def assert_ring(dref, sourcename, filename):
     for i in xrange(li):
         sd.data[0][ind[i]] = 10 * i + 1
     sdref = sref.discretize(topo)
-    ic = topo.mesh.iCompute
+    ic = topo.mesh.compute_index
     assert np.allclose(sd.data[0][ic], sdref.data[0][ic])
 
 
@@ -194,7 +194,7 @@ def assert_ringpole(dref, sourcename, filename):
     for i in xrange(li):
         sd.data[0][ind[i]] = 10 * i + 1
     sdref = sref.discretize(topo)
-    ic = topo.mesh.iCompute
+    ic = topo.mesh.compute_index
     assert np.allclose(sd.data[0][ic], sdref.data[0][ic])
 
 
diff --git a/hysop/domain/tests/test_regular_subset.py b/hysop/domain/tests/test_regular_subset.py
index f20de501d..78c161c9e 100755
--- a/hysop/domain/tests/test_regular_subset.py
+++ b/hysop/domain/tests/test_regular_subset.py
@@ -215,4 +215,4 @@ def test_full2d_in3d():
     i1 = cub.integrate_dfield_allc(vd)
     assert (np.abs(i1 - vref) < 1e-6).all()
     print topo
-    print cub.mesh[topo].iCompute
+    print cub.mesh[topo].compute_index
diff --git a/hysop/domain/tests/test_submesh.py b/hysop/domain/tests/test_submesh.py
index 9f427d7bc..b4efae725 100755
--- a/hysop/domain/tests/test_submesh.py
+++ b/hysop/domain/tests/test_submesh.py
@@ -47,7 +47,7 @@ def check_submesh(discr):
 
         r2 = mesh.convert2local([s for s
                                  in Utils.intersl(gpos, mesh.position)])
-        assert subm.iCompute == r2
+        assert subm.compute_index == r2
         pos = subm.position_in_parent
         subpos = [slice(pos[d].start - gstart[d],
                         pos[d].stop - gstart[d]) for d in xrange(dim)]
diff --git a/hysop/domain/tests/test_subset.py b/hysop/domain/tests/test_subset.py
index 00bb7bcc7..8501fc3b4 100755
--- a/hysop/domain/tests/test_subset.py
+++ b/hysop/domain/tests/test_subset.py
@@ -58,7 +58,7 @@ def check_subset(discr, fileref, class_name):
     scal = Field(domain=topo.domain, is_vector=False, name='s0')
     sd = scal.discretize(topo)
     sd[0][subs.ind[topo][0]] = 2.
-    icompute = topo.mesh.iCompute
+    icompute = topo.mesh.compute_index
     return np.allclose(sd[0][icompute], vd[0][icompute])
 
 
@@ -120,7 +120,7 @@ def init_obst_list(discr, fileref):
     sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
     indices = Subset.union(obs_list, topo)
     sd[0][indices] = 2.
-    icompute = topo.mesh.iCompute
+    icompute = topo.mesh.compute_index
     return np.allclose(sd[0][icompute], vd[0][icompute])
 
 
@@ -145,7 +145,7 @@ def init_subtract(discr, fileref):
     sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
     indices = Subset.subtract(box, s1, topo)
     sd[0][indices] = 2.
-    icompute = topo.mesh.iCompute
+    icompute = topo.mesh.compute_index
     return np.allclose(sd[0][icompute], vd[0][icompute])
 
 
@@ -183,7 +183,7 @@ def init_subtract_lists(discr, fileref):
     sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
     indices = Subset.subtract_list_of_sets(box, obs, topo)
     sd[0][indices] = 2.
-    return np.allclose(sd[0][topo.mesh.iCompute], vd[0][topo.mesh.iCompute])
+    return np.allclose(sd[0][topo.mesh.compute_index], vd[0][topo.mesh.compute_index])
 
 
 def test_subtract_list_3d():
diff --git a/hysop/fields/discrete.py b/hysop/fields/discrete.py
index 803abb2e9..264e1f5e4 100755
--- a/hysop/fields/discrete.py
+++ b/hysop/fields/discrete.py
@@ -192,7 +192,7 @@ class DiscreteField(object):
         """
         result = npw.zeros((self.nb_components))
         gresult = npw.zeros((self.nb_components))
-        ind_ng = self.topology.mesh.iCompute
+        ind_ng = self.topology.mesh.compute_index
         for d in range(self.nb_components):
             result[d] = la.norm(self.data[d][ind_ng]) ** 2
         self.topology.comm.Allreduce(result, gresult)
@@ -207,7 +207,7 @@ class DiscreteField(object):
         """
         result = npw.zeros((self.nb_components))
         gresult = npw.zeros((self.nb_components))
-        ind_ng = self.topology.mesh.iCompute
+        ind_ng = self.topology.mesh.compute_index
         step = self.topology.mesh.space_step
         for d in range(self.nb_components):
             result[d] = (step[d] * la.norm(self.data[d][ind_ng])) ** 2
diff --git a/hysop/gpu/tests/test_gpu_multiresolution_filter.py b/hysop/gpu/tests/test_gpu_multiresolution_filter.py
index b44201bbd..111b4e919 100755
--- a/hysop/gpu/tests/test_gpu_multiresolution_filter.py
+++ b/hysop/gpu/tests/test_gpu_multiresolution_filter.py
@@ -56,11 +56,11 @@ def test_filter_linear():
         f_out.wait()
         valid = [npw.zeros(f_out[0].shape), ]
         valid = func(valid, *topo_coarse.mesh.coords)
-        assert np.allclose(valid[0][topo_coarse.mesh.iCompute],
-                           f_out[0][topo_coarse.mesh.iCompute],
+        assert np.allclose(valid[0][topo_coarse.mesh.compute_index],
+                           f_out[0][topo_coarse.mesh.compute_index],
                            atol=1e-4, rtol=1e-3), \
-            np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                          f_out[0][topo_coarse.mesh.iCompute]))
+            np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                          f_out[0][topo_coarse.mesh.compute_index]))
         if PY_COMPARE:
             f_py = Field(box, formula=func, name='fpy')
             op_py = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
@@ -72,10 +72,10 @@ def test_filter_linear():
             f_py.initialize(topo=topo_fine)
             op_py.apply(simu)
             valid = f_py.discreteFields[topo_coarse]
-            assert np.allclose(valid[0][topo_coarse.mesh.iCompute],
-                               f_out[0][topo_coarse.mesh.iCompute]), \
-                np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                              f_out[0][topo_coarse.mesh.iCompute]))
+            assert np.allclose(valid[0][topo_coarse.mesh.compute_index],
+                               f_out[0][topo_coarse.mesh.compute_index]), \
+                np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                              f_out[0][topo_coarse.mesh.compute_index]))
 
 
 def test_filter_L2_1():
@@ -108,10 +108,10 @@ def test_filter_L2_1():
         f_out.wait()
         valid = [npw.zeros(f_out[0].shape), ]
         valid = func(valid, *topo_coarse.mesh.coords)
-        assert np.allclose(valid[0][topo_coarse.mesh.iCompute],
-                           f_out[0][topo_coarse.mesh.iCompute]), \
-            np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                          f_out[0][topo_coarse.mesh.iCompute]))
+        assert np.allclose(valid[0][topo_coarse.mesh.compute_index],
+                           f_out[0][topo_coarse.mesh.compute_index]), \
+            np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                          f_out[0][topo_coarse.mesh.compute_index]))
         if PY_COMPARE:
             f_py = Field(box, formula=func, name='fpy')
             op_py = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
@@ -122,7 +122,7 @@ def test_filter_L2_1():
             f_py.initialize(topo=topo_fine)
             op_py.apply(simu)
             valid = f_py.discreteFields[topo_coarse]
-            assert np.allclose(valid[0][topo_coarse.mesh.iCompute],
-                               f_out[0][topo_coarse.mesh.iCompute]), \
-                np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                              f_out[0][topo_coarse.mesh.iCompute]))
+            assert np.allclose(valid[0][topo_coarse.mesh.compute_index],
+                               f_out[0][topo_coarse.mesh.compute_index]), \
+                np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                              f_out[0][topo_coarse.mesh.compute_index]))
diff --git a/hysop/mpi/bridge_inter.py b/hysop/mpi/bridge_inter.py
index ecc0af52f..85f37b8a2 100644
--- a/hysop/mpi/bridge_inter.py
+++ b/hysop/mpi/bridge_inter.py
@@ -116,8 +116,8 @@ class BridgeInter(object):
                 self.comm.bcast(current_indices, root=MPI.PROC_NULL)
 
         # Convert numpy arrays to dict of slices ...
-        current_indices = Utils.arrayToDict(current_indices)
-        remote_indices = Utils.arrayToDict(remote_indices)
+        current_indices = Utils.array_to_dict(current_indices)
+        remote_indices = Utils.array_to_dict(remote_indices)
 
         return current_indices, remote_indices
 
diff --git a/hysop/mpi/topology.py b/hysop/mpi/topology.py
index 1ab614715..11a8057b1 100644
--- a/hysop/mpi/topology.py
+++ b/hysop/mpi/topology.py
@@ -508,7 +508,7 @@ class TopoTools(object):
                         root=root)
 
         if toslice:
-            return Utils.arrayToDict(iglob_res)
+            return Utils.array_to_dict(iglob_res)
         else:
             return iglob_res
 
@@ -572,7 +572,7 @@ class TopoTools(object):
                 comm.Gather([iglob[:, rank], MPI.INT], [iglob_res, MPI.INT],
                             root=root)
             if toslice:
-                return Utils.arrayToDict(iglob_res)
+                return Utils.array_to_dict(iglob_res)
             else:
                 return iglob_res
 
diff --git a/hysop/numerics/differential_operations.py b/hysop/numerics/differential_operations.py
index b6096a301..3f78337a6 100755
--- a/hysop/numerics/differential_operations.py
+++ b/hysop/numerics/differential_operations.py
@@ -51,7 +51,7 @@ class DifferentialOperation(object):
         indices : list of slices, optional
             Represents the local mesh on which the operation
             will be applied,
-            like iCompute in :class:`~hysop.domain.mesh.Mesh`.
+            like compute_index in :class:`~hysop.domain.mesh.Mesh`.
             See details in notes.
         reduce_output_shape : boolean, optional
             True to return the result in a reduced array. See notes below.
@@ -68,7 +68,7 @@ class DifferentialOperation(object):
         ** 1 - either invar and outvar are arrays of the same shape, and
         outvar[...] = op(invar[...]) will be computed. In that case
         indices and reduce_output_shape are not required.
-        Indices will be set to topo.mesh.iCompute.
+        Indices will be set to topo.mesh.compute_index.
         ** 2 - or outvar is a smaller array than invar, and
         outvar[...] = op(invar(indices))
         To use case 2, set reduce_output_shape = True and provide indices.
@@ -88,7 +88,7 @@ class DifferentialOperation(object):
 
         self.method = method
         if indices is None:
-            indices = topo.mesh.iCompute
+            indices = topo.mesh.compute_index
             reduce_output_shape = False
         self._indices = indices
         self.fd_scheme = self._init_fd_method(topo, reduce_output_shape)
diff --git a/hysop/numerics/finite_differences.py b/hysop/numerics/finite_differences.py
index c724b4cfc..cecf24579 100644
--- a/hysop/numerics/finite_differences.py
+++ b/hysop/numerics/finite_differences.py
@@ -13,7 +13,7 @@ class FiniteDifference(object):
     Usage :
 
     >> step = topo.mesh_space_step
-    >> scheme = FD_C_4(step, topo.mesh.iCompute)
+    >> scheme = FD_C_4(step, topo.mesh.compute_index)
 
     For a given numpy array (obviously discretized on the topo
     used to compute indices), to compute result as the derivative
@@ -48,7 +48,7 @@ class FiniteDifference(object):
             resolution of the mesh
         indices : list of slices
             Represents the local mesh on which finite-differences
-            will be applied, like iCompute in :class:`~hysop.domain.mesh.Mesh`.
+            will be applied, like compute_index in :class:`~hysop.domain.mesh.Mesh`.
         reduce_output_shape : boolean, optional
             True to return the result in a reduced array. See notes below.
 
diff --git a/hysop/numerics/tests/test_diffOp.py b/hysop/numerics/tests/test_diffOp.py
index 529542e01..555ea48a3 100755
--- a/hysop/numerics/tests/test_diffOp.py
+++ b/hysop/numerics/tests/test_diffOp.py
@@ -108,8 +108,8 @@ def init(discr, vform, wform):
 # Init 2D and 3D
 v3d, w3d, topo3, work3 = init(discr3D, computeVel, computeVort)
 v2d, w2d, topo2, work2 = init(discr2D, computeVel2, computeVort2)
-ic3 = topo3.mesh.iCompute
-ic2 = topo2.mesh.iCompute
+ic3 = topo3.mesh.compute_index
+ic2 = topo2.mesh.compute_index
 
 
 def build_op(op_class, len_result, topo, work):
diff --git a/hysop/numerics/tests/test_update_ghosts.py b/hysop/numerics/tests/test_update_ghosts.py
index 8ec8f843f..b4cfaf2ee 100755
--- a/hysop/numerics/tests/test_update_ghosts.py
+++ b/hysop/numerics/tests/test_update_ghosts.py
@@ -15,7 +15,7 @@ def _setup(ghosts, topo_dim):
     df = f.discreteFields[topo]
     for i, dfd in enumerate(df.data):
         dfd[...] = 0.
-        dfd[topo.mesh.iCompute] = 1. * (i + 1)
+        dfd[topo.mesh.compute_index] = 1. * (i + 1)
     return df, topo
 
 
diff --git a/hysop/numerics/utils.py b/hysop/numerics/utils.py
index 2399a503c..34cd0b8ea 100644
--- a/hysop/numerics/utils.py
+++ b/hysop/numerics/utils.py
@@ -19,7 +19,7 @@ class Utils(object):
         y : list of numpy arrays
             represents a discrete field
         sl : list of slices
-            mesh points indices (like topo.mesh.iCompute)
+            mesh points indices (like topo.mesh.compute_index)
         work: numpy array
             temporary buffer
 
diff --git a/hysop/operator/discrete/baroclinic.py b/hysop/operator/discrete/baroclinic.py
index b02527fb2..607fa5597 100644
--- a/hysop/operator/discrete/baroclinic.py
+++ b/hysop/operator/discrete/baroclinic.py
@@ -58,10 +58,10 @@ class Baroclinic(DiscreteOperator):
 
         self._laplacian = diff_op.Laplacian(
             self.velocity.topology,
-            indices=self.velocity.topology.mesh.iCompute)
+            indices=self.velocity.topology.mesh.compute_index)
         self._gradOp = diff_op.GradS(
             self.velocity.topology,
-            indices=self.velocity.topology.mesh.iCompute,
+            indices=self.velocity.topology.mesh.compute_index,
             method=self.method[SpaceDiscretisation])
 
         # Gravity vector
@@ -73,9 +73,9 @@ class Baroclinic(DiscreteOperator):
     def initialize_velocity(self):
         """Initialize the temporary array 'result' with the velocity"""
         topo = self.velocity.topology
-        iCompute = topo.mesh.iCompute
+        compute_index = topo.mesh.compute_index
         for d in xrange(self.velocity.dimension):
-            self._result[d][iCompute] = -self.velocity[d][iCompute]
+            self._result[d][compute_index] = -self.velocity[d][compute_index]
 
     @debug
     def apply(self, simulation=None):
@@ -95,14 +95,14 @@ class Baroclinic(DiscreteOperator):
         self._synchronizeRho(self.density.data)
 
         topo = self.velocity.topology
-        iCompute = topo.mesh.iCompute
+        compute_index = topo.mesh.compute_index
 
         # result = du/dt = (u^(n)-u^(n-1))/dt
         # result has been initialized with -u^(n-1)
         # and _old_dt equals to the previous dt
         for d in xrange(self.velocity.dimension):
-            self._result[d][iCompute] += self.velocity[d][iCompute]
-            self._result[d][iCompute] /= self._old_dt
+            self._result[d][compute_index] += self.velocity[d][compute_index]
+            self._result[d][compute_index] /= self._old_dt
 
         # result = result + (u . grad)u
         # (u. grad)u = (u.du/dx + v.du/dy + w.du/dz ;
@@ -114,64 +114,64 @@ class Baroclinic(DiscreteOperator):
         # result[X] = result[X] + ((u. grad)u)[X]
         #           = result[X] + u.du/dx + v.du/dy + w.du/dz
         for d in xrange(self.velocity.dimension):
-            self._result[XDIR][iCompute] += \
-                self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+            self._result[XDIR][compute_index] += \
+                self.velocity[d][compute_index] * self._tempGrad[d][compute_index]
 
         self._tempGrad = self._gradOp(
             self.velocity[YDIR:YDIR + 1], self._tempGrad)
         # result[Y] = result[Y] + ((u. grad)u)[Y]
         #           = result[Y] + u.dv/dx + v.dv/dy + w.dv/dz
         for d in xrange(self.velocity.dimension):
-            self._result[YDIR][iCompute] += \
-                self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+            self._result[YDIR][compute_index] += \
+                self.velocity[d][compute_index] * self._tempGrad[d][compute_index]
 
         self._tempGrad = self._gradOp(
             self.velocity[ZDIR:ZDIR + 1], self._tempGrad)
         # result[Z] = result[Z] + ((u. grad)u)[Z]
         #           = result[Z] + u.dw/dx + v.dw/dy + w.dw/dz
         for d in xrange(self.velocity.dimension):
-            self._result[ZDIR][iCompute] += \
-                self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+            self._result[ZDIR][compute_index] += \
+                self.velocity[d][compute_index] * self._tempGrad[d][compute_index]
 
         # result = result - nu*\Laplacian u (-g) = gradP/rho
         for d in xrange(self.velocity.dimension):
             self._tempGrad[d:d + 1] = self._laplacian(
                 self.velocity[d:d + 1], self._tempGrad[d:d + 1])
         for d in xrange(self.velocity.dimension):
-            self._tempGrad[d][iCompute] *= self.viscosity
-            self._result[d][iCompute] -= self._tempGrad[d][iCompute]
+            self._tempGrad[d][compute_index] *= self.viscosity
+            self._result[d][compute_index] -= self._tempGrad[d][compute_index]
 
         # gravity term : result = result - g
         for d in xrange(self.velocity.dimension):
-            self._result[2][iCompute] -= self._gravity[d]
+            self._result[2][compute_index] -= self._gravity[d]
 
         # baroclinicTerm = -(gradRho/rho) x (gradP/rho)
         self._tempGrad = self._gradOp(self.density[0:1], self._tempGrad)
         for d in xrange(self.velocity.dimension):
-            self._tempGrad[d][iCompute] = \
-                self._tempGrad[d][iCompute] / self.density[0][iCompute]
-
-        self._baroclinicTerm[0][iCompute] = \
-            - self._tempGrad[1][iCompute] * self._result[2][iCompute]
-        self._baroclinicTerm[0][iCompute] += \
-            self._tempGrad[2][iCompute] * self._result[1][iCompute]
-        self._baroclinicTerm[1][iCompute] = \
-            - self._tempGrad[2][iCompute] * self._result[0][iCompute]
-        self._baroclinicTerm[1][iCompute] += \
-            self._tempGrad[0][iCompute] * self._result[2][iCompute]
-        self._baroclinicTerm[2][iCompute] = \
-            - self._tempGrad[0][iCompute] * self._result[1][iCompute]
-        self._baroclinicTerm[2][iCompute] += \
-            self._tempGrad[1][iCompute] * self._result[0][iCompute]
+            self._tempGrad[d][compute_index] = \
+                self._tempGrad[d][compute_index] / self.density[0][compute_index]
+
+        self._baroclinicTerm[0][compute_index] = \
+            - self._tempGrad[1][compute_index] * self._result[2][compute_index]
+        self._baroclinicTerm[0][compute_index] += \
+            self._tempGrad[2][compute_index] * self._result[1][compute_index]
+        self._baroclinicTerm[1][compute_index] = \
+            - self._tempGrad[2][compute_index] * self._result[0][compute_index]
+        self._baroclinicTerm[1][compute_index] += \
+            self._tempGrad[0][compute_index] * self._result[2][compute_index]
+        self._baroclinicTerm[2][compute_index] = \
+            - self._tempGrad[0][compute_index] * self._result[1][compute_index]
+        self._baroclinicTerm[2][compute_index] += \
+            self._tempGrad[1][compute_index] * self._result[0][compute_index]
 
         # vorti(n+1) = vorti(n) + dt * baroclinicTerm
         for d in xrange(self.vorticity.dimension):
-            self._baroclinicTerm[d][iCompute] *= dt
-            self.vorticity[d][iCompute] += self._baroclinicTerm[d][iCompute]
+            self._baroclinicTerm[d][compute_index] *= dt
+            self.vorticity[d][compute_index] += self._baroclinicTerm[d][compute_index]
 
 
         # reinitialise for next iteration
         # velo(n-1) update
         for d in xrange(self.velocity.dimension):
-            self._result[d][iCompute] = -self.velocity.data[d][iCompute]
+            self._result[d][compute_index] = -self.velocity.data[d][compute_index]
         self._old_dt = dt
diff --git a/hysop/operator/discrete/baroclinic_from_rhs.py b/hysop/operator/discrete/baroclinic_from_rhs.py
index c503edbd6..1434833bc 100644
--- a/hysop/operator/discrete/baroclinic_from_rhs.py
+++ b/hysop/operator/discrete/baroclinic_from_rhs.py
@@ -54,8 +54,8 @@ class BaroclinicFromRHS(DiscreteOperator):
         #self._synchronizeRHS(self.rhs.data)
 
         topo = self.vorticity.topology
-        iCompute = topo.mesh.iCompute
+        compute_index = topo.mesh.compute_index
 
         # vorti(n+1) = vorti(n) + dt * baroclinicTerm
         for d in xrange(self.vorticity.dimension):
-            self.vorticity[d][iCompute] += self.rhs[d][iCompute] * dt
+            self.vorticity[d][compute_index] += self.rhs[d][compute_index] * dt
diff --git a/hysop/operator/discrete/density.py b/hysop/operator/discrete/density.py
index ffea3a72b..d72df5ca7 100644
--- a/hysop/operator/discrete/density.py
+++ b/hysop/operator/discrete/density.py
@@ -42,14 +42,14 @@ class DensityVisco_d(DiscreteOperator):
         assert simulation is not None, \
             "Missing simulation value for computation."
 
-        iCompute = self.density.topology.mesh.iCompute
+        compute_index = self.density.topology.mesh.compute_index
 
         # Density reconstruction
-        if self.density[0][iCompute].all() <= np.absolute(
+        if self.density[0][compute_index].all() <= np.absolute(
                 self.densityVal[1] - self.densityVal[0]) / 2.0:
-            self.density[0][iCompute] = self.densityVal[1]
+            self.density[0][compute_index] = self.densityVal[1]
         else:
-            self.density[0][iCompute] = self.densityVal[0]
+            self.density[0][compute_index] = self.densityVal[0]
 
         # Viscosity reconstruction :
         # nu = nu1 + (nu2 - nu1) * (density - rho1)/(rho2 - rho1)
diff --git a/hysop/operator/discrete/energy_enstrophy.py b/hysop/operator/discrete/energy_enstrophy.py
index b5101215f..e22512fca 100644
--- a/hysop/operator/discrete/energy_enstrophy.py
+++ b/hysop/operator/discrete/energy_enstrophy.py
@@ -50,8 +50,8 @@ class EnergyEnstrophy(DiscreteOperator):
 
     def _set_work_arrays(self, rwork=None, iwork=None):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        w_ind = self.vorticity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_v = self.velocity.data[0][v_ind].shape
         shape_w = self.velocity.data[0][w_ind].shape
         # setup for rwork, iwork is useless.
@@ -75,8 +75,8 @@ class EnergyEnstrophy(DiscreteOperator):
 
     def get_work_properties(self):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        w_ind = self.vorticity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_v = self.velocity.data[0][v_ind].shape
         shape_w = self.velocity.data[0][w_ind].shape
         if shape_v == shape_w:
@@ -94,7 +94,7 @@ class EnergyEnstrophy(DiscreteOperator):
         vd = self.velocity
         # get the list of computation points (no ghosts)
         nbc = vd.nb_components
-        v_ind = self.velocity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
         # Integrate (locally) velocity ** 2
         local_energy = 0.
         for i in xrange(nbc):
@@ -104,7 +104,7 @@ class EnergyEnstrophy(DiscreteOperator):
         # --- Enstrophy computation ---
         vortd = self.vorticity
         nbc = vortd.nb_components
-        w_ind = self.vorticity.topology.mesh.iCompute
+        w_ind = self.vorticity.topology.mesh.compute_index
         # Integrate (locally) vorticity ** 2
         work = self._rwork[-1]
         local_enstrophy = 0.
diff --git a/hysop/operator/discrete/monitoringPoints.py b/hysop/operator/discrete/monitoringPoints.py
index 74738987b..de8ae0e8b 100644
--- a/hysop/operator/discrete/monitoringPoints.py
+++ b/hysop/operator/discrete/monitoringPoints.py
@@ -35,7 +35,7 @@ class MonitoringPoints(DiscreteOperator):
         super(MonitoringPoints, self).__init__(variables=[velocity, vorticity],
                                        **kwds)
         topo_v = self.velocity.topology
-        self.shape_v = self.velocity.data[0][topo_v.mesh.iCompute].shape
+        self.shape_v = self.velocity.data[0][topo_v.mesh.compute_index].shape
         self.space_step = topo_v.mesh.space_step
         self.length = topo_v.domain.length
         self.origin = topo_v.domain.origin
@@ -53,8 +53,8 @@ class MonitoringPoints(DiscreteOperator):
 
     def _set_work_arrays(self, rwork=None, iwork=None):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        w_ind = self.vorticity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_v = self.velocity.data[0][v_ind].shape
         shape_w = self.velocity.data[0][w_ind].shape
         # setup for rwork, iwork is useless.
@@ -78,8 +78,8 @@ class MonitoringPoints(DiscreteOperator):
 
     def get_work_properties(self):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        w_ind = self.vorticity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_v = self.velocity.data[0][v_ind].shape
         shape_w = self.velocity.data[0][w_ind].shape
         if shape_v == shape_w:
diff --git a/hysop/operator/discrete/multiphase_gradp.py b/hysop/operator/discrete/multiphase_gradp.py
index 0a2697b62..b15e780f5 100644
--- a/hysop/operator/discrete/multiphase_gradp.py
+++ b/hysop/operator/discrete/multiphase_gradp.py
@@ -54,10 +54,10 @@ class GradP(DiscreteOperator):
 
         self._laplacian = diff_op.Laplacian(
             self.velocity.topology,
-            indices=self.velocity.topology.mesh.iCompute)
+            indices=self.velocity.topology.mesh.compute_index)
         self._gradOp = diff_op.GradS(
             self.velocity.topology,
-            indices=self.velocity.topology.mesh.iCompute,
+            indices=self.velocity.topology.mesh.compute_index,
             method=self.method[SpaceDiscretisation])
 
         # Gravity vector
@@ -70,9 +70,9 @@ class GradP(DiscreteOperator):
     def initialize_velocity(self):
         """Initialize the temporary array 'result' with the velocity"""
         topo = self.velocity.topology
-        iCompute = topo.mesh.iCompute
+        compute_index = topo.mesh.compute_index
         for d in xrange(self.velocity.dimension):
-            self._result[d][iCompute] = -self.velocity[d][iCompute]
+            self._result[d][compute_index] = -self.velocity[d][compute_index]
 
     @debug
     def apply(self, simulation=None):
@@ -89,14 +89,14 @@ class GradP(DiscreteOperator):
         self._synchronizeVel(self.velocity.data)
 
         topo = self.velocity.topology
-        iCompute = topo.mesh.iCompute
+        compute_index = topo.mesh.compute_index
 
         # result = du/dt = (u^(n)-u^(n-1))/dt
         # result has been initialized with -u^(n-1)
         # and _old_dt equals to the previous dt
         for d in xrange(self.velocity.dimension):
-            self._result[d][iCompute] += self.velocity[d][iCompute]
-            self._result[d][iCompute] /= self._old_dt
+            self._result[d][compute_index] += self.velocity[d][compute_index]
+            self._result[d][compute_index] /= self._old_dt
 
         # result = result + (u . grad)u
         # (u. grad)u = (u.du/dx + v.du/dy + w.du/dz ;
@@ -108,36 +108,36 @@ class GradP(DiscreteOperator):
         # result[X] = result[X] + ((u. grad)u)[X]
         #           = result[X] + u.du/dx + v.du/dy + w.du/dz
         for d in xrange(self.velocity.dimension):
-            self._result[XDIR][iCompute] += \
-                self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+            self._result[XDIR][compute_index] += \
+                self.velocity[d][compute_index] * self._tempGrad[d][compute_index]
 
         self._tempGrad = self._gradOp(
             self.velocity[YDIR:YDIR + 1], self._tempGrad)
         # result[Y] = result[Y] + ((u. grad)u)[Y]
         #           = result[Y] + u.dv/dx + v.dv/dy + w.dv/dz
         for d in xrange(self.velocity.dimension):
-            self._result[YDIR][iCompute] += \
-                self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+            self._result[YDIR][compute_index] += \
+                self.velocity[d][compute_index] * self._tempGrad[d][compute_index]
 
         self._tempGrad = self._gradOp(
             self.velocity[ZDIR:ZDIR + 1], self._tempGrad)
         # result[Z] = result[Z] + ((u. grad)u)[Z]
         #           = result[Z] + u.dw/dx + v.dw/dy + w.dw/dz
         for d in xrange(self.velocity.dimension):
-            self._result[ZDIR][iCompute] += \
-                self.velocity[d][iCompute] * self._tempGrad[d][iCompute]
+            self._result[ZDIR][compute_index] += \
+                self.velocity[d][compute_index] * self._tempGrad[d][compute_index]
 
         # result = result - nu*\Laplacian u (-g) = gradP/rho
         for d in xrange(self.velocity.dimension):
             self._tempGrad[d:d + 1] = self._laplacian(
                 self.velocity[d:d + 1], self._tempGrad[d:d + 1])
         for d in xrange(self.velocity.dimension):
-            self._tempGrad[d][iCompute] *= self.viscosity
-            self._result[d][iCompute] -= self._tempGrad[d][iCompute]
+            self._tempGrad[d][compute_index] *= self.viscosity
+            self._result[d][compute_index] -= self._tempGrad[d][compute_index]
 
         # gravity term : result = result - g
         for d in xrange(self.velocity.dimension):
-            self._result[2][iCompute] -= self._gravity[d]
+            self._result[2][compute_index] -= self._gravity[d]
 
         for d in xrange(self.velocity.dimension):
             self.gradp[d][...] = self._result[d]
@@ -145,5 +145,5 @@ class GradP(DiscreteOperator):
         # reinitialise for next iteration
         # velo(n-1) update
         for d in xrange(self.velocity.dimension):
-            self._result[d][iCompute] = -self.velocity.data[d][iCompute]
+            self._result[d][compute_index] = -self.velocity.data[d][compute_index]
         self._old_dt = dt
diff --git a/hysop/operator/discrete/multiresolution_filter.py b/hysop/operator/discrete/multiresolution_filter.py
index b35f3229c..4639aa25b 100644
--- a/hysop/operator/discrete/multiresolution_filter.py
+++ b/hysop/operator/discrete/multiresolution_filter.py
@@ -113,19 +113,19 @@ class FilterFineToCoarse(DiscreteOperator):
             for i_y in xrange(len(self._rmsh_method.weights)):
                 for i_z in xrange(len(self._rmsh_method.weights)):
                     self._sl_coarse.append((
-                        slice(self._mesh_out.iCompute[0].start -
+                        slice(self._mesh_out.compute_index[0].start -
                               self._rmsh_method.shift + i_x,
-                              self._mesh_out.iCompute[0].stop -
+                              self._mesh_out.compute_index[0].stop -
                               self._rmsh_method.shift + i_x,
                               None),
-                        slice(self._mesh_out.iCompute[1].start -
+                        slice(self._mesh_out.compute_index[1].start -
                               self._rmsh_method.shift + i_y,
-                              self._mesh_out.iCompute[1].stop -
+                              self._mesh_out.compute_index[1].stop -
                               self._rmsh_method.shift + i_y,
                               None),
-                        slice(self._mesh_out.iCompute[2].start -
+                        slice(self._mesh_out.compute_index[2].start -
                               self._rmsh_method.shift + i_z,
-                              self._mesh_out.iCompute[2].stop -
+                              self._mesh_out.compute_index[2].stop -
                               self._rmsh_method.shift + i_z,
                               None)
                     ))
diff --git a/hysop/operator/discrete/profiles.py b/hysop/operator/discrete/profiles.py
index 10c175f3d..74ca03e4f 100644
--- a/hysop/operator/discrete/profiles.py
+++ b/hysop/operator/discrete/profiles.py
@@ -41,7 +41,7 @@ class Profiles(DiscreteOperator):
         super(Profiles, self).__init__(variables=[velocity, vorticity],
                                        **kwds)
         topo_v = self.velocity.topology
-        self.shape_v = self.velocity.data[0][topo_v.mesh.iCompute].shape
+        self.shape_v = self.velocity.data[0][topo_v.mesh.compute_index].shape
         self.space_step = topo_v.mesh.space_step
         self.length = topo_v.domain.length
         self.origin = topo_v.domain.origin
@@ -64,8 +64,8 @@ class Profiles(DiscreteOperator):
 
     def _set_work_arrays(self, rwork=None, iwork=None):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        w_ind = self.vorticity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_v = self.velocity.data[0][v_ind].shape
         shape_w = self.velocity.data[0][w_ind].shape
         # setup for rwork, iwork is useless.
@@ -89,8 +89,8 @@ class Profiles(DiscreteOperator):
 
     def get_work_properties(self):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        w_ind = self.vorticity.topology.mesh.iCompute
+        v_ind = self.velocity.topology.mesh.compute_index
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_v = self.velocity.data[0][v_ind].shape
         shape_w = self.velocity.data[0][w_ind].shape
         if shape_v == shape_w:
diff --git a/hysop/operator/discrete/residual.py b/hysop/operator/discrete/residual.py
index 498924fd0..2fa568b12 100644
--- a/hysop/operator/discrete/residual.py
+++ b/hysop/operator/discrete/residual.py
@@ -26,7 +26,7 @@ class Residual(DiscreteOperator):
         assert 'variables' not in kwds, 'variables parameter is useless.'
         super(Residual, self).__init__(variables=[vorticity], **kwds)
         topo_w = self.vorticity.topology
-        self.shape_w = self.vorticity.data[0][topo_w.mesh.iCompute].shape
+        self.shape_w = self.vorticity.data[0][topo_w.mesh.compute_index].shape
         self.space_step = topo_w.mesh.space_step
         self.length = topo_w.domain.length
         self.origin = topo_w.domain.origin
@@ -41,7 +41,7 @@ class Residual(DiscreteOperator):
 
     def _set_work_arrays(self, rwork=None, iwork=None):
 
-        w_ind = self.vorticity.topology.mesh.iCompute
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_w = self.vorticity.data[0][w_ind].shape
         # setup for rwork, iwork is useless.
         if rwork is None:
@@ -56,13 +56,13 @@ class Residual(DiscreteOperator):
 
     def get_work_properties(self):
 
-        w_ind = self.vorticity.topology.mesh.iCompute
+        w_ind = self.vorticity.topology.mesh.compute_index
         shape_w = self.vorticity.data[0][w_ind].shape
         return {'rwork': [shape_w], 'iwork': None}
 
     def initialize_vortPrev(self):
         
-        w_ind = self.vorticity.topology.mesh.iCompute
+        w_ind = self.vorticity.topology.mesh.compute_index
         for d in xrange(self.vorticity.dimension):
             self._vortPrev[d][w_ind] = self.vorticity[d][w_ind]
 
@@ -85,7 +85,7 @@ class Residual(DiscreteOperator):
         # get the list of computation points (no ghosts)
         wd = self.vorticity
         nbc = wd.nb_components
-        w_ind = self.vorticity.topology.mesh.iCompute
+        w_ind = self.vorticity.topology.mesh.compute_index
         for i in xrange(nbc):
             self._rwork[0][...] = (wd[i][w_ind] -
                                    self._vortPrev[i][w_ind]) ** 2
diff --git a/hysop/operator/discrete/velocity_correction.py b/hysop/operator/discrete/velocity_correction.py
index a720d6960..0053ecd9b 100755
--- a/hysop/operator/discrete/velocity_correction.py
+++ b/hysop/operator/discrete/velocity_correction.py
@@ -131,10 +131,10 @@ class VelocityCorrection_D(DiscreteOperator):
         # Computation of the required velocity shift
         # for the current state
         self.computeCorrection()
-        iCompute = self.topo.mesh.iCompute
+        compute_index = self.topo.mesh.compute_index
 
         # Apply shift to velocity
-        self.velocity[XDIR][iCompute] += self.velocity_shift[XDIR]
+        self.velocity[XDIR][compute_index] += self.velocity_shift[XDIR]
         start = self.velocity.nb_components
         # reminder : self.rates =[vx_shift, flowrates[Y], flowrate[Z],
         #                         vort_mean[X], vort_mean[Y], vort_mean[Z]]
diff --git a/hysop/operator/energy_enstrophy.py b/hysop/operator/energy_enstrophy.py
index 4f0bd403a..a59c1edeb 100644
--- a/hysop/operator/energy_enstrophy.py
+++ b/hysop/operator/energy_enstrophy.py
@@ -52,8 +52,8 @@ class EnergyEnstrophy(Computational):
             raise RuntimeError(msg)
         vd = self.discreteFields[self.velocity]
         wd = self.discreteFields[self.vorticity]
-        v_ind = vd.topology.mesh.iCompute
-        w_ind = wd.topology.mesh.iCompute
+        v_ind = vd.topology.mesh.compute_index
+        w_ind = wd.topology.mesh.compute_index
         shape_v = vd[0][v_ind].shape
         shape_w = wd[0][w_ind].shape
         if shape_v == shape_w:
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index e0771b289..5653b0aa3 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -123,7 +123,7 @@ class HDF_IO(Computational):
         # Global resolution for hdf5 output (warning : this must
         # be the whole domain resolution, not the subset resolution)
         self._global_resolution = list(refmesh.global_resolution())
-        self._slices = refmesh.iCompute
+        self._slices = refmesh.compute_index
         if refmesh.on_proc:
             sl = list(refmesh.position)
         else:
diff --git a/hysop/operator/monitoringPoints.py b/hysop/operator/monitoringPoints.py
index 3066fda28..d67adc99f 100644
--- a/hysop/operator/monitoringPoints.py
+++ b/hysop/operator/monitoringPoints.py
@@ -45,8 +45,8 @@ class MonitoringPoints(Computational):
             raise RuntimeError(msg)
         vd = self.discreteFields[self.velocity]
         wd = self.discreteFields[self.vorticity]
-        v_ind = vd.topology.mesh.iCompute
-        w_ind = wd.topology.mesh.iCompute
+        v_ind = vd.topology.mesh.compute_index
+        w_ind = wd.topology.mesh.compute_index
         shape_v = vd[0][v_ind].shape
         shape_w = wd[0][w_ind].shape
         if shape_v == shape_w:
diff --git a/hysop/operator/profiles.py b/hysop/operator/profiles.py
index a7ffb9e30..ac7dfbdfe 100644
--- a/hysop/operator/profiles.py
+++ b/hysop/operator/profiles.py
@@ -50,8 +50,8 @@ class Profiles(Computational):
             raise RuntimeError(msg)
         vd = self.discreteFields[self.velocity]
         wd = self.discreteFields[self.vorticity]
-        v_ind = vd.topology.mesh.iCompute
-        w_ind = wd.topology.mesh.iCompute
+        v_ind = vd.topology.mesh.compute_index
+        w_ind = wd.topology.mesh.compute_index
         shape_v = vd[0][v_ind].shape
         shape_w = wd[0][w_ind].shape
         if shape_v == shape_w:
diff --git a/hysop/operator/residual.py b/hysop/operator/residual.py
index 1a5de67ae..077f52632 100644
--- a/hysop/operator/residual.py
+++ b/hysop/operator/residual.py
@@ -34,7 +34,7 @@ class Residual(Computational):
             msg += 'before any call to this function.'
             raise RuntimeError(msg)
         wd = self.discreteFields[self.vorticity]
-        w_ind = wd.topology.mesh.iCompute
+        w_ind = wd.topology.mesh.compute_index
         shape_w = wd[0][w_ind].shape
         return {'rwork': [shape_w], 'iwork': None}
 
diff --git a/hysop/operator/tests/test_differential.py b/hysop/operator/tests/test_differential.py
index 7a217638b..2ab50a59a 100755
--- a/hysop/operator/tests/test_differential.py
+++ b/hysop/operator/tests/test_differential.py
@@ -103,7 +103,7 @@ def check(op, ref_formula, topo, op_dim=3, order=4):
     res_d = result.discreteFields[topo]
 
     # Compare results with reference
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
     err = topo.mesh.space_step ** order
     dim = topo.domain.dimension
     for i in xrange(result.nb_components):
diff --git a/hysop/operator/tests/test_hdf5_io.py b/hysop/operator/tests/test_hdf5_io.py
index 368836fd2..3e338fbae 100755
--- a/hysop/operator/tests/test_hdf5_io.py
+++ b/hysop/operator/tests/test_hdf5_io.py
@@ -98,7 +98,7 @@ def test_write_read_scalar_3D():
     reader.setup()
     sc3d = scal3D.discretize(topo)
     scref = scalRef.discretize(topo)
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
     for d in xrange(scal3D.nb_components):
         sc3d.data[d][...] = 0.0
         assert not np.allclose(scref.data[d][ind], sc3d.data[d][ind])
@@ -138,7 +138,7 @@ def test_write_read_scalar_3D_defaults():
 
     sc3d = scal3D.discretize(topo)
     scref = scalRef.discretize(topo)
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
     for d in xrange(scal3D.nb_components):
         sc3d.data[d][...] = scref.data[d][...]
         scref.data[d][...] = 0.0
@@ -194,7 +194,7 @@ def test_write_read_vectors_3D_defaults():
 
     v3d = velo.discretize(topo)
     w3d = vorti.discretize(topo)
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
 
     # Copy current values of v3 and w3 into buff1 and buff2, for comparison
     # after reader.apply, below.
@@ -266,7 +266,7 @@ def test_write_read_vectors_3D():
 
     v3d = velo.discretize(topo)
     w3d = vorti.discretize(topo)
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
 
     buff1 = Field(domain=dom, name='buff1', is_vector=True)
     buff2 = Field(domain=dom, name='buff2', is_vector=True)
@@ -323,8 +323,8 @@ def test_write_read_subset_1():
     assert os.path.exists(fullpath + '_00001.h5')
 
     v3d = velo.discretize(topo)
-    ind = topo.mesh.iCompute
-    indsubset = mybox.mesh[topo].iCompute
+    ind = topo.mesh.compute_index
+    indsubset = mybox.mesh[topo].compute_index
 
     buff1 = Field(domain=dom, name='buff1', is_vector=True)
 
@@ -378,8 +378,8 @@ def test_write_read_subset_2():
     assert os.path.exists(fullpath + '_00001.h5')
 
     v3d = velo.discretize(topo)
-    ind = topo.mesh.iCompute
-    indsubset = mybox.mesh[topo].iCompute
+    ind = topo.mesh.compute_index
+    indsubset = mybox.mesh[topo].compute_index
 
     buff1 = Field(domain=dom, name='buff1', is_vector=True)
 
diff --git a/hysop/operator/tests/test_multiphase_gradp.py b/hysop/operator/tests/test_multiphase_gradp.py
index 477c9e48c..f9a3a000c 100755
--- a/hysop/operator/tests/test_multiphase_gradp.py
+++ b/hysop/operator/tests/test_multiphase_gradp.py
@@ -59,7 +59,7 @@ def test_gradp():
     op.discretize()
     op.setup()
     topo = op.discreteFields[velo].topology
-    iCompute = topo.mesh.iCompute
+    compute_index = topo.mesh.compute_index
     velo.initialize(topo=topo)
     op.initialize_velocity()
     res.initialize(topo=topo)
@@ -67,6 +67,6 @@ def test_gradp():
     op.apply(simu)
     d_res = res.discreteFields[topo]
     d_true_res = true_res.discreteFields[topo]
-    assert np.allclose(d_res[0][iCompute], d_true_res[0][iCompute], atol=8e-3)
-    assert np.allclose(d_res[1][iCompute], d_true_res[1][iCompute], atol=8e-3)
-    assert np.allclose(d_res[2][iCompute], d_true_res[2][iCompute])
+    assert np.allclose(d_res[0][compute_index], d_true_res[0][compute_index], atol=8e-3)
+    assert np.allclose(d_res[1][compute_index], d_true_res[1][compute_index], atol=8e-3)
+    assert np.allclose(d_res[2][compute_index], d_true_res[2][compute_index])
diff --git a/hysop/operator/tests/test_multiresolution_filter.py b/hysop/operator/tests/test_multiresolution_filter.py
index 44e3e5e7e..fbda089f3 100755
--- a/hysop/operator/tests/test_multiresolution_filter.py
+++ b/hysop/operator/tests/test_multiresolution_filter.py
@@ -55,11 +55,11 @@ def filter(d_fine, d_coarse, func, method, atol=1e-8, rtol=1e-5):
     op.apply(simu)
     valid = [npw.zeros(f_out[0].shape), ]
     valid = func(valid, *topo_coarse.mesh.coords)
-    e = np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                      f_out[0][topo_coarse.mesh.iCompute]))
-    err = atol + rtol * np.max(np.abs(valid[0][topo_coarse.mesh.iCompute]))
-    return np.allclose(f_out[0][topo_coarse.mesh.iCompute],
-                       valid[0][topo_coarse.mesh.iCompute],
+    e = np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                      f_out[0][topo_coarse.mesh.compute_index]))
+    err = atol + rtol * np.max(np.abs(valid[0][topo_coarse.mesh.compute_index]))
+    return np.allclose(f_out[0][topo_coarse.mesh.compute_index],
+                       valid[0][topo_coarse.mesh.compute_index],
                        atol=atol, rtol=rtol), e, err
 
 
@@ -188,11 +188,11 @@ def test_filter_linear():
     op.apply(simu)
     valid = [npw.zeros(f_out[0].shape), ]
     valid = func(valid, *topo_coarse.mesh.coords)
-    assert np.allclose(valid[0][topo_coarse.mesh.iCompute],
-                       f_out[0][topo_coarse.mesh.iCompute],
+    assert np.allclose(valid[0][topo_coarse.mesh.compute_index],
+                       f_out[0][topo_coarse.mesh.compute_index],
                        atol=1e-4, rtol=1e-3), \
-        np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                      f_out[0][topo_coarse.mesh.iCompute]))
+        np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                      f_out[0][topo_coarse.mesh.compute_index]))
 
 
 def test_filter_l2_1():
@@ -214,8 +214,8 @@ def test_filter_l2_1():
     op.apply(simu)
     valid = [npw.zeros(f_out[0].shape), ]
     valid = func(valid, *topo_coarse.mesh.coords)
-    assert np.allclose(valid[0][topo_coarse.mesh.iCompute],
-                       f_out[0][topo_coarse.mesh.iCompute]), \
-        np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
-                      f_out[0][topo_coarse.mesh.iCompute]))
+    assert np.allclose(valid[0][topo_coarse.mesh.compute_index],
+                       f_out[0][topo_coarse.mesh.compute_index]), \
+        np.max(np.abs(valid[0][topo_coarse.mesh.compute_index] -
+                      f_out[0][topo_coarse.mesh.compute_index]))
 
diff --git a/hysop/operator/tests/test_penalization.py b/hysop/operator/tests/test_penalization.py
index 9026e986b..75bd69718 100755
--- a/hysop/operator/tests/test_penalization.py
+++ b/hysop/operator/tests/test_penalization.py
@@ -100,7 +100,7 @@ def check_penal(penal, sref, vref, scal, velo):
     sd = scal.discretize(topo)
     simu = Simulation(nbIter=3)
     penal.apply(simu)
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
     assert np.allclose(sd.data[0][ind], sref.data[0][ind])
     for d in xrange(vd.nb_components):
         assert np.allclose(vd.data[d][ind], vref.data[d][ind])
@@ -282,7 +282,7 @@ def check_penal_vort(penal, wref, vref, vorti, velo):
     velo.initialize(topo=topo)
     vd = velo.discretize(topo)
     wd = vorti.discretize(topo)
-    ind = topo.mesh.iCompute
+    ind = topo.mesh.compute_index
 
     simu = Simulation(nbIter=200)
     penal.apply(simu)
diff --git a/hysop/operator/tests/test_redistribute.py b/hysop/operator/tests/test_redistribute.py
index a8bae9037..f5450b9a1 100755
--- a/hysop/operator/tests/test_redistribute.py
+++ b/hysop/operator/tests/test_redistribute.py
@@ -76,7 +76,7 @@ def test_distribute_intra_1():
     red.apply()
     red.wait()
     vd_2 = velocity.discretize(default_topo)
-    ind_ng = default_topo.mesh.iCompute
+    ind_ng = default_topo.mesh.compute_index
     for d in xrange(velocity.nb_components):
         assert np.allclose(vd_2.data[d][ind_ng], wd.data[d][ind_ng])
 
@@ -136,7 +136,7 @@ def test_distribute_intra_2():
     red.apply()
     red.wait()
     vd_2 = velocity.discretize(default_topo)
-    ind_ng = default_topo.mesh.iCompute
+    ind_ng = default_topo.mesh.compute_index
     wd = vorticity.discretize(default_topo)
     for d in xrange(velocity.nb_components):
         assert np.allclose(vd_2.data[d][ind_ng], wd.data[d][ind_ng])
@@ -182,7 +182,7 @@ def test_distribute_intra_3():
     red.apply()
     red.wait()
     vd_2 = velocity.discretize(target_topo)
-    ind_ng = target_topo.mesh.iCompute
+    ind_ng = target_topo.mesh.compute_index
 
     for d in xrange(velocity.nb_components):
         assert np.allclose(vd_2.data[d][ind_ng], wd.data[d][ind_ng])
@@ -231,7 +231,7 @@ def test_distribute_intra_4():
     red.apply()
     red.wait()
     vd_2 = velocity.discretize(target_topo)
-    ind_ng = target_topo.mesh.iCompute
+    ind_ng = target_topo.mesh.compute_index
 
     for d in xrange(velocity.nb_components):
         assert np.allclose(vd_2.data[d][ind_ng], wd.data[d][ind_ng])
@@ -284,7 +284,7 @@ def test_distribute_intra_5():
     red.wait()
     vd = velocity.discretize(target_topo)
     wd = vorticity.discretize(target_topo)
-    ind_ng = target_topo.mesh.iCompute
+    ind_ng = target_topo.mesh.compute_index
 
     for d in xrange(velocity.nb_components):
         assert np.allclose(vd.data[d][ind_ng], rd.data[d][ind_ng])
@@ -385,7 +385,7 @@ def test_distribute_overlap():
     red.wait()
     if target_topo is not None:
         vd_2 = velocity.discretize(target_topo)
-        ind_ng = target_topo.mesh.iCompute
+        ind_ng = target_topo.mesh.compute_index
         for d in xrange(velocity.nb_components):
             assert np.allclose(vd_2.data[d][ind_ng], wd.data[d][ind_ng])
 
@@ -483,7 +483,7 @@ def test_distribute_inter_2():
         vd = velocity.discretize(topo=topo_CPU)
         wd = reference.discretize(topo=topo_CPU)
         vnorm = velocity.norm(topo_CPU)
-        ind = topo_CPU.mesh.iCompute
+        ind = topo_CPU.mesh.compute_index
         wnorm = reference.norm(topo_CPU)
         assert np.allclose(vnorm, wnorm)
         for d in xrange(dom.dimension):
@@ -527,7 +527,7 @@ def test_distribute_inter_3():
         wd = reference.discretize(topo=topo2)
         vd = velocity.discretize(topo=topo2)
         vnorm = velocity.norm(topo2)
-        ind = topo2.mesh.iCompute
+        ind = topo2.mesh.compute_index
         assert np.allclose(vnorm, wnorm)
         for d in xrange(dom.dimension):
             assert np.allclose(wd.data[d][ind], vd.data[d][ind])
@@ -567,7 +567,7 @@ def test_distribute_inter_4():
     if dom_tasks.is_on_task(CPU):
         wd = reference.discretize(topo=topo)
         vd = velocity.discretize(topo=topo)
-        ind = topo.mesh.iCompute
+        ind = topo.mesh.compute_index
         for d in xrange(dom.dimension):
             assert np.allclose(wd.data[d][ind], vd.data[d][ind])
 
-- 
GitLab