diff --git a/examples/NSDebug_faux2D.py b/examples/NSDebug_faux2D.py
index a6540feed2bd024932eb177dee1e6164fb08778f..16902a65893e27ebbf3bca7c9b3a28332f0cfb23 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 7b0b93cc9da4cdf96cf906913b1d9f6fad09b226..0a08866abc350b39154f6411dea6a20914b08b81 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 e2c7504a45fb3adc744f6cdc29da4a6ebd784b87..a7b10a5e153ee138eed02fab6b7abeaafc6e5e0d 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 4b0ec75762da0b9b3d3cb7e18bdeb4256c1f8615..01ea25d8677820689dbbc0a7d40c8ab176c4f171 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 46b0735fe14b66e1663be28b64948eae743640ee..412f5fb63175e4ab7e23f872405f90c23215ea3c 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 02a36beb2437b34ec9f2c057e4d203d75c4ee5e3..5dbdd9b245907ebe8768901806919a9c88d6aa02 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 0ae2b4868b213a6358f65b93656c441287e3b0dd..261abb444423a07bbd0ac43020e35a782a3f8d1a 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 3ae7c79c9d72dc47bc3d661176c46ee86523b03f..39dd869dd276e04531a3ef3814bc2c1ac651f066 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 f20de501d21ede246456ed100600a1b7ca40857a..78c161c9ebc781d06f8fb1c59a009c2459a5b3e4 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 9f427d7bc3c45af14cc1ccdfab13f623ea17db09..b4efae725daad40d11b7d4cc6c51e570ae6a9070 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 00bb7bcc7e9ec03c4f2888aa2ca2d739b6cbad15..8501fc3b41600e3b9a4efb28eae4338d08ae069e 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 803abb2e9f6acd4d75f67812c1b7d7b6aa219287..264e1f5e41f0a79c4b2fba8fdf5f4d2275abe1ab 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 b44201bbd0bf6b5969e445303196cdf58a5dd827..111b4e919d50a313e37a4d56779a6b51af3a51f8 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 ecc0af52fac0d9fc7b6edc4b4f758692e3f10027..85f37b8a289b4ebb4eaa45c720206f3e0fbc8027 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 1ab614715ff7bb37cd467ca003727506349d634c..11a8057b16f61a46a24bfbbe901358e65f381486 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 b6096a30118c23b9872d7bb1865eefa74e28622c..3f78337a60ee3fdc12309aa4c4406d0f26fc759c 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 c724b4cfc8c7b913384e20a18f06bb16ec70c1d8..cecf245794feca8d17fd730f4395b4a4bf5cb341 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 529542e012b024bd7dd0f7144a96affa79c4f156..555ea48a3df2931084495ef8b4d421c3ec268224 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 8ec8f843fff29cbd6227d5e6c128463d4f8cc747..b4cfaf2ee219b5729c2d26d594a25338f51f99a2 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 2399a503ca15496f8dd9563a7b2366ae0a4860be..34cd0b8eaff68a6beb79d76b8dbf64fd21ba75dd 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 b02527fb21f4b62ea9f2947b8490e587f42638f8..607fa5597081f84bf3c6dcdf7345a8ff91f38ebc 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 c503edbd68807e3ee5beb200add79f47efc474c3..1434833bcef1be7d4c3ea4476a7e1cfca2d98cac 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 ffea3a72b31bc01482405a04d7dc37327f5dcf85..d72df5ca7e9e7ca7d98f99c2d4dc8402af3704ff 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 b5101215f80fcc17b1aa98d04ce19900b4d64492..e22512fcaa54cd9861948cf6199284de83e52d3a 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 74738987b7cebf540b8c6a6423f473795fa8eb55..de8ae0e8bff4c13c8eaa0b4a3361003f43041073 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 0a2697b623b9e5f79b1147dbffb6c96e92fb654f..b15e780f5bc1e829d1bc38ffcdc711a413b47ed1 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 b35f3229c4520b25fede521e191ec3c287475981..4639aa25b23190a2f26ac06136e84952ccbb3b87 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 10c175f3df043fafd6e028fa1be5f0098126d22a..74ca03e4f0157b8d4f421650b7a5908ca1aba115 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 498924fd0b4f63be944f6dd385c0a4cdebd156fb..2fa568b12ae65f6fceff5d0210f4d9f0f41ef41a 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 a720d6960a4f12cba5dadad59309b6467642c945..0053ecd9b8542868ef3884d7cfdefa89c8f3ef2a 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 4f0bd403adf9cfcd87d2dfdcec8dcf654ba2cac8..a59c1edeb1816023ceb6113400be8ffc0eb0a3e3 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 e0771b2890c6098baebd3e0ca1572d75643d9b05..5653b0aa322a27b6ccd36b1f8ffc2da5cffff974 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 3066fda287f34011417866f37a5095c649fe3fc1..d67adc99fc418be061980f985aca2f1561d02f7c 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 a7ffb9e30b2027dabcb45159f76cce133921cb76..ac7dfbdfe8486f5a6b0695f48dade759752a87b6 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 1a5de67ae04acdec7e08051b6b36fc8aeb601fb0..077f52632b343d85a6f096231bcf7aac7332115e 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 7a217638b23f2335873563e391bdc67dc301f9ff..2ab50a59af60b251ff00bf1df7ff428f866b7aab 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 368836fd23275fc59a5f0d2178d0a8953001f071..3e338fbaeec4d3da38d3f79ee4e62d51daaf98dd 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 477c9e48c28fe66acf895d77eb9b86bd1ffbfa98..f9a3a000ce7f1f35c703263b89bb2412234ffc0f 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 44e3e5e7ec5c22522d6ad7cbc98c72640a48209f..fbda089f3a0ff9acd524a30aa6e312532527079c 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 9026e986bf6eacb419f4b4f11798ab343fed7d99..75bd6971806ba7d1097437680f66069646f9260e 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 a8bae9037aa321bed99017c317f03e9cfc706cdc..f5450b9a1d3d5a9f2eb12a0b2cac8aed855c0edd 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])