Repository: PMEAL/OpenPNM
Pull request #430 (2015-08-21 22:10:47)

Subdivide modification
The bug  #428 is fixed now. It was a minor bug due to the wrong naming inside the method. 
Also some topological methods and their unit tests have been added.
I had these changes in one of my private branches for a long time. There are still some changes for improving subdivide method which are important, (especially for multi-scale systems), but since I do not have the necessary tests for them at this moment, I'll make another pull request for them in near future.


Patch:
diff --git a/OpenPNM/Network/__Cubic__.py b/OpenPNM/Network/__Cubic__.py
index d194fd984..8bfed418d 100644
--- a/OpenPNM/Network/__Cubic__.py
+++ b/OpenPNM/Network/__Cubic__.py
@@ -337,3 +337,46 @@ class Cubic(GenericNetwork):
     def subdivide(self, pores=[], shape=[], labels=[]):
         topo.subdivide(network=self, pores=pores, shape=shape, labels=labels)
     subdivide.__doc__ = topo.subdivide.__doc__
+
+    def template_sphere_shell(outer_radius=None, inner_radius=0):
+        r"""
+        This method generates an image array of a sphere shell for a cubic network.
+
+        Parameters
+        ----------
+        outer_radius : array_like
+        Number of the nodes in the outer radius of the shell
+
+        inner_radius : float
+        Number of the nodes in the inner radius of the shell
+
+        """
+
+        if outer_radius is None:
+            raise Exception('No outer radius has been sent!')
+        if inner_radius is None:
+            raise Exception('Number of nodes in the inner radius cannot be '
+                            'None!')
+        rmax = np.array(outer_radius, ndmin=1)
+        rmin = np.array(inner_radius, ndmin=1)
+        s_rmax = np.size(rmax)
+        s_rmin = np.size(rmin)
+        if not ((s_rmax in [1, 3]) and (s_rmin in [1, 3])):
+            raise Exception('In this method, each radius can be scalar or '
+                            'array with components along all xyz directions.')
+        s_u_rmax = np.size(np.unique(rmax))
+        s_u_rmin = np.size(np.unique(rmin))
+        if not ((s_u_rmax == 1) and (s_u_rmin == 1)):
+            raise Exception('In this method, all components of radius should '
+                            'be unique values along all xyz directions.')
+        pnum = 2 * np.ones(3) * rmax - 1
+        Rx, Ry, Rz = np.array(pnum, dtype=np.int32)
+        x, y, z = np.indices((Rx, Ry, Rz))
+        x = x - (Rx - 1)/2
+        y = y - (Ry - 1)/2
+        z = z - (Rz - 1)/2
+        img = x ** 2 + y ** 2 + z ** 2 < np.unique(rmax) ** 2
+        if not np.all(rmin == 0):
+            img_min = x ** 2 + y ** 2 + z ** 2 > np.unique(rmin) ** 2
+            img = img * img_min
+        return (img)
diff --git a/OpenPNM/Utilities/__topology__.py b/OpenPNM/Utilities/__topology__.py
index ea8a17152..1226ab01c 100644
--- a/OpenPNM/Utilities/__topology__.py
+++ b/OpenPNM/Utilities/__topology__.py
@@ -6,6 +6,7 @@ Network.tools.topology: Assorted topological manipulation methods
 
 """
 import scipy as _sp
+import numpy as _np
 import scipy.sparse as _sprs
 import scipy.spatial as _sptl
 from OpenPNM.Base import logging as _logging
@@ -351,7 +352,6 @@ class topology(object):
         if method == 'nearest':
             P1 = P_network
             P2 = P_donor + N_init['pore']  # Increment pores on donor
-            P = _sp.hstack((P1, P2))
             C1 = network['pore.coords'][P_network]
             C2 = donor['pore.coords'][P_donor]
             D = _sp.spatial.distance.cdist(C1, C2)
@@ -412,8 +412,8 @@ class topology(object):
         Notes
         -----
         It creates the connections in a format which is acceptable by
-        the default OpenPNM connection key ('throat.conns') and adds them to the
-        network.
+        the default OpenPNM connection key ('throat.conns') and adds them to
+        the network.
 
         Examples
         --------
@@ -440,6 +440,28 @@ class topology(object):
         conns = _sp.vstack([array1, array2]).T
         self.extend(network=network, throat_conns=conns, labels=labels)
 
+    def find_centroid(coords=None):
+        r'''
+        It finds the coordinates of the centroid of the sent pores.
+        '''
+        l = _np.float64(len(coords))
+        x, y, z = coords.T
+        sx = _np.sum(x)
+        sy = _np.sum(y)
+        sz = _np.sum(z)
+        c = _np.array([sx/l, sy/l, sz/l], ndmin=1)
+        return c
+
+    def find_pores_distance(network, pores1=None, pores2=None):
+        r'''
+        It finds the distance between two group of pores.
+        '''
+        from scipy.spatial.distance import cdist
+        p1 = _sp.array(pores1, ndmin=1)
+        p2 = _sp.array(pores2, ndmin=1)
+        coords = network['pore.coords']
+        return cdist(coords[p1], coords[p2])
+
     def subdivide(self, network, pores, shape, labels=[]):
         r'''
         It trim the pores and replace them by cubic networks with the sent shape.
@@ -525,8 +547,8 @@ class topology(object):
             else:
                 for ind in [0, 1]:
                     loc = (non_single_labels[ind] == l)
-                    pores = new_net.pores(non_single_labels[ind][loc])
-                    new_net['pore.surface_' + l][pores] = True
+                    temp_pores = new_net.pores(non_single_labels[ind][loc])
+                    new_net['pore.surface_' + l][temp_pores] = True
 
         old_coords = _sp.copy(new_net['pore.coords'])
         if labels == []:


Test patch:
diff --git a/test/integration/test_topological_manipulations.py b/test/integration/test_topological_manipulations.py
index 40e6a433a..2d6a2346c 100644
--- a/test/integration/test_topological_manipulations.py
+++ b/test/integration/test_topological_manipulations.py
@@ -15,8 +15,6 @@ def test_subdivide():
     pn.subdivide(pores=nano_pores, shape=[4, 4, 4], labels='nano')
     assert pn.Np == (125+4*64-4)
     assert pn.Nt == (300+(4*144)-16+15*16+16)
-    ctrl.export(network=pn, filename=join(TEMP_DIR, 'nano'))
-
 
 def test_clone_and_trim():
     ctrl.clear()
@@ -66,3 +64,33 @@ def test_stitch():
     assert pn.Np == 3*pn2.Np  # Ensure number of pores increased again
     assert pn.Nt == (3*pn2.Nt + 2*Nx*Ny)  # Ensure correct num of new throats
     ctrl.clear()
+
+
+def test_distance_center():
+    shape = sp.array([7, 5, 9])
+    spacing = sp.array([2, 1, 0.5])
+    pn = OpenPNM.Network.Cubic(shape=shape, spacing=spacing)
+    sx, sy, sz = spacing
+    center_coord = sp.around(topology.find_centroid(pn['pore.coords']), 7)
+    cx, cy, cz = center_coord
+    coords = pn['pore.coords']
+    x, y, z = coords.T
+    coords = sp.concatenate((coords, center_coord.reshape((1, 3))))
+    pn['pore.center'] = False
+    mask1 = (x <= (cx + sx/2)) * (y <= (cy + sy/2)) * (z <= (cz + sz/2))
+    mask2 = (x >= (cx - sx/2)) * (y >= (cy - sy/2)) * (z >= (cz - sz/2))
+    center_pores_mask = pn.Ps[mask1 * mask2]
+    pn['pore.center'][center_pores_mask] = True
+    center = pn.Ps[pn['pore.center']]
+    L1 = sp.amax(topology.find_pores_distance(network=pn,
+                                              pores1=center,
+                                              pores2=pn.Ps))
+    L2 = sp.amax(topology.find_pores_distance(network=pn,
+                                              pores1=pn.Ps,
+                                              pores2=pn.Ps))
+    l1 = ((shape[0] - 1) * sx) ** 2
+    l2 = ((shape[1] - 1) * sy) ** 2
+    l3 = ((shape[2] - 1) * sz) ** 2
+    L3 = sp.sqrt(l1 + l2 + l3)
+    assert sp.around(L1 * 2, 7) == sp.around(L2, 7)
+    assert sp.around(L2, 7) == sp.around(L3, 7)
diff --git a/test/unit/Network/CubicTest.py b/test/unit/Network/CubicTest.py
index 9b06bd3da..c7d863ddb 100644
--- a/test/unit/Network/CubicTest.py
+++ b/test/unit/Network/CubicTest.py
@@ -92,3 +92,19 @@ class CubicTest:
         L = self.net.domain_length(face_1=self.net.pores('top'),
                                    face_2=self.net.pores('bottom'))
         assert sp.allclose(L, 4, rtol=1e-02)
+
+    def test_template_sphere_shell(self):
+        from OpenPNM.Utilities import topology
+        spacing = sp.array([0.5])
+        r_o = [5, 5, 5]
+        r_in = [2, 2, 2]
+        img = OpenPNM.Network.Cubic.template_sphere_shell(outer_radius=r_o,
+                                                          inner_radius=r_in)
+        pn_sphere = OpenPNM.Network.Cubic(template=img, spacing=spacing)
+        assert pn_sphere.Np == 452
+        img2 = OpenPNM.Network.Cubic.template_sphere_shell(outer_radius=r_o)
+        pn_sphere = OpenPNM.Network.Cubic(template=img2, spacing=spacing)
+        L1 = sp.amax(topology.find_pores_distance(network=pn_sphere,
+                                                  pores1=pn_sphere.Ps,
+                                                  pores2=pn_sphere.Ps))
+        assert L1 < sp.sqrt(3) * 8 * 0.5