diff --git a/src/kanapy/initializations.py b/src/kanapy/initializations.py
index 0e281692..7e602c0a 100644
--- a/src/kanapy/initializations.py
+++ b/src/kanapy/initializations.py
@@ -503,7 +503,7 @@ def create_voxels(self, sim_box):
         node_count = 0
         elmt_count = 0
         # loop over the duplicate pairs
-        for (mk, nk), (mj, nj), (mi, ni) in itertools.product(pointsX_dup, pointsY_dup, pointsZ_dup):
+        for (mi, ni), (mj, nj), (mk, nk) in itertools.product(pointsX_dup, pointsY_dup, pointsZ_dup):
 
             # Find the center of each voxel and update the center dictionary
             elmt_count += 1
diff --git a/src/kanapy/input_output.py b/src/kanapy/input_output.py
index 731cdc7b..22947695 100644
--- a/src/kanapy/input_output.py
+++ b/src/kanapy/input_output.py
@@ -424,7 +424,7 @@ def import_voxels(file, path='./'):
     phases = np.zeros(nvox, dtype=int)
     grain_dict = dict()
     grain_phase_dict = dict()
-    gr_arr = grains.flatten(order='F')
+    gr_arr = grains.flatten(order='C')
     if 'Grains' in data.keys():
         if 'Orientation' in data['Grains'][grain_keys[-1]].keys():
             grain_ori_dict = dict()
@@ -489,7 +489,7 @@ def import_voxels(file, path='./'):
     ms.mesh.grains = grains
     ms.mesh.grain_dict = grain_dict
     ms.mesh.grain_ori_dict = grain_ori_dict
-    ms.mesh.phases = phases.reshape(sh, order='F')
+    ms.mesh.phases = phases.reshape(sh, order='C')
     ms.mesh.grain_phase_dict = grain_phase_dict
     ms.mesh.ngrains_phase = ngrain
     if 0 in ms.mesh.grain_dict.keys():
diff --git a/src/kanapy/plotting.py b/src/kanapy/plotting.py
index 1534d9ea..66194c7e 100644
--- a/src/kanapy/plotting.py
+++ b/src/kanapy/plotting.py
@@ -61,6 +61,8 @@ def plot_voxels_3D(data, Ngr=1, sliced=False, dual_phase=False,
         vox_b = mask
     if dual_phase:
         Ngr = 2
+    else:
+        Ngr = np.max(data)
     if clist is None:
         cm = plt.cm.get_cmap(cmap, Ngr)
         colors = cm(data.astype(int))
diff --git a/src/kanapy/textures.py b/src/kanapy/textures.py
index 9d7def69..d7c5b746 100644
--- a/src/kanapy/textures.py
+++ b/src/kanapy/textures.py
@@ -152,14 +152,14 @@ def __init__(self, fname, matname=None, gs_min=3, vf_min=0.03, plot=True):
                     eng.hold('on', nargout=0)
                     eng.mtexColorbar
                 except:
-                    logging.warning('ODF too large for plotting')'''
+                    logging.warning('ODF too large for plotting')
                 # workaround:
                 eng.workspace["ori"] = data['ori']
                 eng.workspace["cs"] = data['cs']
                 try:
                     eng.eval(f"plotPDF(ori,Miller(0,0,1,cs),'points','all')", nargout=0)
                 except:
-                    logging.warning('Plotting of orientations failed.')
+                    logging.warning('Plotting of orientations failed.')'''
             # Evaluate grain shape statistics
             # generate dict for statistical input for geometry module
 
diff --git a/src/kanapy/voxelization.py b/src/kanapy/voxelization.py
index b796c6a9..a283948d 100644
--- a/src/kanapy/voxelization.py
+++ b/src/kanapy/voxelization.py
@@ -50,8 +50,8 @@ def assign_voxels_to_ellipsoid(cooDict, Ellipsoids, voxel_dict, vf_target=None):
     Nvox = len(test_ids)
     test_points = np.array(list(cooDict.values()))
 
-    # array defining ellipsoid growth for each stage of while loop
-    growth = iter(list(np.arange(1.0, 15.0, 0.1)))
+    # scale factor for defining ellipsoid growth for each stage of while loop
+    scale = 1.0
     remaining_voxels = set(list(cooDict.keys()))
     assigned_voxels = set()
     vf_cur = 0.
@@ -63,7 +63,7 @@ def assign_voxels_to_ellipsoid(cooDict, Ellipsoids, voxel_dict, vf_target=None):
 
     while vf_cur < vf_target:
         # call the growth value for the ellipsoids
-        scale = next(growth)
+        scale += 0.1
         for ellipsoid in Ellipsoids:
             # scale ellipsoid dimensions by the growth factor and generate surface points
             ellipsoid.a, ellipsoid.b, ellipsoid.c = scale * ellipsoid.a, scale * ellipsoid.b, scale * ellipsoid.c
@@ -379,7 +379,7 @@ def poly2vox():
     for igr, vlist in mesh.grain_dict.items():
         vlist = np.array(vlist) - 1
         gr_arr[vlist] = igr
-    mesh.grains = np.reshape(gr_arr, mesh.dim, order='F')
+    mesh.grains = np.reshape(gr_arr, mesh.dim, order='C')
 
     # generate array of voxelized structure with phase numbers
     # and dict of phase numbers for each grain
@@ -395,7 +395,7 @@ def poly2vox():
         mesh.ngrains_phase[ip] += 1
     ind = np.nonzero(ph_arr < 0.0)[0]
     ph_arr[ind] = 1  # assign phase 1 to empty voxels
-    mesh.phases = np.reshape(ph_arr, mesh.dim, order='F')
+    mesh.phases = np.reshape(ph_arr, mesh.dim, order='C')
     vf_cur = 1.0 - len(ind) / mesh.nvox
 
     print('Completed RVE voxelization')
diff --git a/tests/test_voxelization.py b/tests/test_voxelization.py
index 35ea2648..4759d288 100644
--- a/tests/test_voxelization.py
+++ b/tests/test_voxelization.py
@@ -213,7 +213,7 @@ def test_voxelizationRoutine():
     mesh.create_voxels(sim_box)
     Particles = particle_generator([parDict], sim_box, None)
     mesh = voxelizationRoutine(Particles, mesh, nphases)
-    assert mesh.vox_center_dict[5][0] == 0.9
-    assert mesh.voxel_dict[7] == [29, 30, 26, 25, 31, 32, 28, 27]
+    assert mesh.vox_center_dict[5][2] == 0.9
+    assert np.all(np.sort(mesh.voxel_dict[7]) == [25, 26, 27, 28, 29, 30, 31, 32])
     assert mesh.grains.shape == (15, 15, 15)