Repository: PMEAL/OpenPNM
Pull request #1156 (2018-12-12 02:44:42)

Fix 1145 find surface 2d


Patch:
diff --git a/openpnm/network/Cubic.py b/openpnm/network/Cubic.py
index 272c4dfba..2265abec7 100644
--- a/openpnm/network/Cubic.py
+++ b/openpnm/network/Cubic.py
@@ -149,27 +149,13 @@ class Cubic(GenericNetwork):
 
         self['pore.coords'] = points
         self['throat.conns'] = pairs
-        # Label faces
-        x, y, z = self['pore.coords'].T
-        labels = ['front', 'back', 'left', 'right', 'bottom', 'top']
-        for label in labels:
-            if 'pore.'+label not in self.keys():
-                self['pore.'+label] = False
-        self['pore.front'][x <= x.min()] = True
-        self['pore.back'][x >= x.max()] = True
-        self['pore.left'][y <= y.min()] = True
-        self['pore.right'][y >= y.max()] = True
-        self['pore.bottom'][z <= z.min()] = True
-        self['pore.top'][z >= z.max()] = True
-        # Label surface pores
-        self['pore.surface'] = False
-        Ps = self.pores(labels)
-        self['pore.surface'][Ps] = True
-        self['throat.surface'] = False
-        Ts = self.find_neighbor_throats(pores=Ps, mode='xnor')
-        self['throat.surface'][Ts] = True
         self['pore.internal'] = True
         self['throat.internal'] = True
+        topotools.label_faces(self)
+        Ps = self.pores('surface')
+        Ts = self.find_neighbor_throats(pores=Ps, mode='xnor')
+        self['throat.surface'] = False
+        self['throat.surface'][Ts] = True
         # Scale network to requested spacing
         self['pore.coords'] *= spacing
 
diff --git a/openpnm/network/CubicTemplate.py b/openpnm/network/CubicTemplate.py
index 11a139e1c..1f5d53d72 100644
--- a/openpnm/network/CubicTemplate.py
+++ b/openpnm/network/CubicTemplate.py
@@ -76,3 +76,8 @@ class CubicTemplate(Cubic):
         self['pore.drop'] = template.flatten() == 0
         topotools.trim(network=self, pores=self.pores('drop'))
         del self['pore.drop']
+        # remove labels pertaining to surface pores, then redo post-trim
+        self.clear(mode='labels')
+        self['pore.internal'] = True
+        self['throat.internal'] = True
+        topotools.find_surface_pores(self)
diff --git a/openpnm/topotools/__init__.py b/openpnm/topotools/__init__.py
index a34a2fad2..d81cfe440 100644
--- a/openpnm/topotools/__init__.py
+++ b/openpnm/topotools/__init__.py
@@ -14,6 +14,7 @@ from .topotools import add_boundary_pores
 from .topotools import bond_percolation
 from .topotools import clone_pores
 from .topotools import connect_pores
+from .topotools import dimensionality
 from .topotools import extend
 from .topotools import find_path
 from .topotools import find_surface_pores
diff --git a/openpnm/topotools/topotools.py b/openpnm/topotools/topotools.py
index 961c559b7..01237b405 100644
--- a/openpnm/topotools/topotools.py
+++ b/openpnm/topotools/topotools.py
@@ -1030,10 +1030,12 @@ def label_faces(network, tol=0.0, label='surface'):
 
     label : string
         An identifying label to isolate the pores on the faces of the network.
-        default is 'surface'.
+        The default is 'surface'.  Surface pores can be found using
+        ``find_surface_pores``.
 
     """
-    if label not in network.labels():
+    label = label.split('.', 1)[-1]
+    if 'pore.'+label not in network.labels():
         find_surface_pores(network, label=label)
     Psurf = network['pore.'+label]
     crds = network['pore.coords']
@@ -1043,12 +1045,16 @@ def label_faces(network, tol=0.0, label='surface'):
     yspan = ymax - ymin
     zmin, zmax = sp.amin(crds[:, 2]), sp.amax(crds[:, 2])
     zspan = zmax - zmin
-    network['pore.back'] = (crds[:, 0] >= (xmax - tol*xspan)) * Psurf
-    network['pore.right'] = (crds[:, 1] >= (ymax - tol*yspan)) * Psurf
-    network['pore.top'] = (crds[:, 2] >= (zmax - tol*zspan)) * Psurf
-    network['pore.front'] = (crds[:, 0] <= (xmin + tol*xspan)) * Psurf
-    network['pore.left'] = (crds[:, 1] <= (ymin + tol*yspan)) * Psurf
-    network['pore.bottom'] = (crds[:, 2] <= (zmin + tol*zspan)) * Psurf
+    dims = dimensionality(network)
+    if dims[0]:
+        network['pore.front'] = (crds[:, 0] <= (xmin + tol*xspan)) * Psurf
+        network['pore.back'] = (crds[:, 0] >= (xmax - tol*xspan)) * Psurf
+    if dims[1]:
+        network['pore.left'] = (crds[:, 1] <= (ymin + tol*yspan)) * Psurf
+        network['pore.right'] = (crds[:, 1] >= (ymax - tol*yspan)) * Psurf
+    if dims[2]:
+        network['pore.top'] = (crds[:, 2] >= (zmax - tol*zspan)) * Psurf
+        network['pore.bottom'] = (crds[:, 2] <= (zmin + tol*zspan)) * Psurf
 
 
 def find_surface_pores(network, markers=None, label='surface'):
@@ -1099,26 +1105,39 @@ def find_surface_pores(network, markers=None, label='surface'):
     >>> net.num_pores(['top','bottom', 'left', 'right', 'front','back'])
     98
 
-    This function is mostly useful for unique networks such as spheres, random
-    topology, or networks that have been subdivied.
-
     """
     import scipy.spatial as sptl
     if markers is None:
-        (xmax, ymax, zmax) = sp.amax(network['pore.coords'], axis=0)
-        (xmin, ymin, zmin) = sp.amin(network['pore.coords'], axis=0)
-        xave = (xmin+xmax)/2
-        yave = (ymin+ymax)/2
-        zave = (zmin+zmax)/2
-        markers = [[xmax + xave, yave, zave],
-                   [xmin - xave, yave, zave],
-                   [xave, ymax + yave, zave],
-                   [xave, ymin - yave, zave],
-                   [xave, yave, zmax + zave],
-                   [xave, yave, zmin - zave]]
-    markers = sp.atleast_2d(markers)
-    tri = sptl.Delaunay(network['pore.coords'], incremental=True)
-    tri.add_points(markers)
+        dims = dimensionality(network)
+        coords = network['pore.coords'][:, dims]
+        # normalize coords to a 1 unit cube centered on origin
+        coords -= sp.amin(coords, axis=0)
+        coords /= sp.amax(coords, axis=0)
+        coords -= 0.5
+        npts = max((network.Np/10, 100))
+        if sum(dims) == 1:
+            network['pore.'+label] = True
+            return
+        if sum(dims) == 2:
+            r = 0.75
+            theta = sp.linspace(0, 2*sp.pi, npts, dtype=float)
+            x = r*sp.cos(theta)
+            y = r*sp.sin(theta)
+            markers = sp.vstack((x, y)).T
+        if sum(dims) == 3:
+            r = 1.00
+            indices = sp.arange(0, npts, dtype=float) + 0.5
+            phi = sp.arccos(1 - 2*indices/npts)
+            theta = sp.pi * (1 + 5**0.5) * indices
+            x = r*sp.cos(theta) * sp.sin(phi)
+            y = r*sp.sin(theta) * sp.sin(phi)
+            z = r*sp.cos(phi)
+            markers = sp.vstack((x, y, z)).T
+    else:
+        coords = network['pore.coords']
+        markers = sp.atleast_2d(markers)
+    pts = sp.vstack((coords, markers))
+    tri = sptl.Delaunay(pts, incremental=False)
     (indices, indptr) = tri.vertex_neighbor_vertices
     for k in range(network.Np, tri.npoints):
         neighbors = indptr[indices[k]:indices[k+1]]
@@ -1129,8 +1148,30 @@ def find_surface_pores(network, markers=None, label='surface'):
         network['pore.'+label][neighbors] = True
 
 
+def dimensionality(network):
+    r"""
+    Checks the dimensionality of the network
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network whose dimensionality is to be checked
+
+    Returns
+    -------
+    Returns an 3-by-1 array containing ``True`` for each axis that contains
+    multiple values, indicating that the pores are spatially distributed
+    in that direction.
+    """
+    v = [True, True, True]
+    for i in [0, 1, 2]:
+        v[i] = not sp.allclose(network['pore.coords'][:, i],
+                               network['pore.coords'][1, i])
+    return v
+
+
 def clone_pores(network, pores, labels=['clone'], mode='parents'):
-    r'''
+    r"""
     Clones the specified pores and adds them to the network
 
     Parameters
@@ -1151,7 +1192,7 @@ def clone_pores(network, pores, labels=['clone'], mode='parents'):
         - 'siblings': Clones are only connected to each other in the same
                       manner as parents were connected
         - 'isolated': No connections between parents or siblings
-    '''
+    """
     if len(network.project.geometries()) > 0:
         logger.warning('Network has active Geometries, new pores must be \
                         assigned a Geometry')
@@ -1470,7 +1511,7 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
 
 def find_pore_to_pore_distance(network, pores1=None, pores2=None):
     r'''
-    Find the distance between all pores on set one to each pore in set 2
+    Find the distance between all pores on set 1 to each pore in set 2
 
     Parameters
     ----------
@@ -1490,6 +1531,13 @@ def find_pore_to_pore_distance(network, pores1=None, pores2=None):
     The distance between pore *i* in ``pores1`` and *j* in ``pores2`` is
     located at *(i, j)* and *(j, i)* in the distance matrix.
 
+    Notes
+    -----
+    This function computes and returns a distance matrix, which is a dense
+    matrix of size Np_1 by Np_2, so can get large.  For distances between
+    larger sets a KD-tree approach would be better, which is available in
+    ``scipy.spatial``.
+
     '''
     from scipy.spatial.distance import cdist
     p1 = sp.array(pores1, ndmin=1)
@@ -1709,33 +1757,33 @@ def merge_pores(network, pores, labels=['merged']):
         len(pores[0])
     except (TypeError, IndexError):
         pores = [pores]
-        
+
     N = len(pores)
     NBs, XYZs = [], []
-    
+
     for Ps in pores:
         NBs.append(network.find_neighbor_pores(pores=Ps,
                                                mode='union',
                                                flatten=True,
                                                include_input=False))
         XYZs.append(network['pore.coords'][Ps].mean(axis=0))
-    
+
     extend(network, pore_coords=XYZs, labels=labels)
     Pnew = network.Ps[-N::]
-    
+
     # Possible throats between new pores: This only happens when running in
     # batch mode, i.e. multiple groups of pores are to be merged. In case
-    # some of these groups share elements, possible throats between the 
+    # some of these groups share elements, possible throats between the
     # intersecting elements is not captured and must be added manually.
     pores_set = [set(items) for items in pores]
     NBs_set = [set(items) for items in NBs]
-    ps1, ps2 = [], []    
+    ps1, ps2 = [], []
     from itertools import combinations
     for i, j in combinations(range(N), 2):
         if not NBs_set[i].isdisjoint(pores_set[j]):
             ps1.append([network.Ps[-N+i]])
             ps2.append([network.Ps[-N+j]])
-    
+
     # Add (possible) connections between the new pores
     connect_pores(network, pores1=ps1, pores2=ps2, labels=labels)
     # Add connections between the new pores and the rest of the network
@@ -1788,7 +1836,7 @@ def _template_sphere_disc(dim, outer_radius, inner_radius):
     return img
 
 
-def template_sphere_shell(outer_radius, inner_radius=0):
+def template_sphere_shell(outer_radius, inner_radius=0, dim=3):
     r"""
     This method generates an image array of a sphere-shell. It is useful for
     passing to Cubic networks as a ``template`` to make spherical shaped
@@ -1803,13 +1851,17 @@ def template_sphere_shell(outer_radius, inner_radius=0):
         Number of nodes in the inner radius of the shell.  a value of 0 will
         result in a solid sphere.
 
+    dim : scalar
+        Controls the number of dimensions of the result.  3 returns a sphere,
+        while 2 returns a disk.
+
     Returns
     -------
     A Numpy array containing 1's to demarcate the sphere-shell, and 0's
     elsewhere.
 
     """
-    img = _template_sphere_disc(dim=3, outer_radius=outer_radius,
+    img = _template_sphere_disc(dim=dim, outer_radius=outer_radius,
                                 inner_radius=inner_radius)
     return img
 


Test patch:
diff --git a/tests/unit/network/CubicDualTest.py b/tests/unit/network/CubicDualTest.py
index 6e78da3b9..d20f766f6 100644
--- a/tests/unit/network/CubicDualTest.py
+++ b/tests/unit/network/CubicDualTest.py
@@ -46,8 +46,6 @@ class CubicDualTest:
         assert net.num_pores('right') == 9
         assert net.num_pores('front') == 9
         assert net.num_pores('back') == 9
-        assert net.num_pores('top') == 41
-        assert net.num_pores('bottom') == 41
         assert net.num_throats('interconnect') == 96
 
     def test_generation_2D_XZ(self):
@@ -55,8 +53,6 @@ class CubicDualTest:
                                    label_2='secondary')
         assert net.Np == 57
         assert net.Nt == 176
-        assert net.num_pores('left') == 41
-        assert net.num_pores('right') == 41
         assert net.num_pores('front') == 9
         assert net.num_pores('back') == 9
         assert net.num_pores('top') == 9
@@ -70,8 +66,6 @@ class CubicDualTest:
         assert net.Nt == 176
         assert net.num_pores('left') == 9
         assert net.num_pores('right') == 9
-        assert net.num_pores('front') == 41
-        assert net.num_pores('back') == 41
         assert net.num_pores('top') == 9
         assert net.num_pores('bottom') == 9
         assert net.num_throats('interconnect') == 96
@@ -85,8 +79,6 @@ class CubicDualTest:
         assert net.num_pores('right') == 9
         assert net.num_pores('front') == 9
         assert net.num_pores('back') == 9
-        assert net.num_pores('top') == 41
-        assert net.num_pores('bottom') == 41
         assert net.num_throats('interconnect') == 96
 
 
diff --git a/tests/unit/network/LabelTest.py b/tests/unit/network/LabelTest.py
index 309834ef4..bf7d5ad87 100644
--- a/tests/unit/network/LabelTest.py
+++ b/tests/unit/network/LabelTest.py
@@ -30,6 +30,19 @@ class LabelTest:
         assert net.num_pores('left') == 16
         assert net.num_pores(['left', 'top'], mode='xnor') == 4
 
+    def test_cubic_2D(self):
+        net = op.network.Cubic(shape=[5, 5, 1])
+        assert 'pore.top' not in net.labels()
+        assert 'pore.bottom' not in net.labels()
+        assert net.num_pores('surface') == 16
+        net = op.network.Cubic(shape=[5, 1, 5])
+        assert 'pore.left' not in net.labels()
+        assert 'pore.right' not in net.labels()
+        assert net.num_pores('surface') == 16
+        net = op.network.Cubic(shape=[1, 5, 5])
+        assert 'pore.front' not in net.labels()
+        assert 'pore.back' not in net.labels()
+        assert net.num_pores('surface') == 16
 
 if __name__ == '__main__':
 
diff --git a/tests/unit/topotools/TopotoolsTest.py b/tests/unit/topotools/TopotoolsTest.py
index f0479890f..4c9c12ca0 100644
--- a/tests/unit/topotools/TopotoolsTest.py
+++ b/tests/unit/topotools/TopotoolsTest.py
@@ -56,7 +56,7 @@ class TopotoolsTest:
         net.clear(mode='labels')
         assert net.labels() == ['pore.all', 'throat.all']
         topotools.find_surface_pores(network=net)
-        assert net.num_pores('surface') == 42
+        assert net.num_pores('surface') == 66
 
     def test_find_pore_to_pore_distance(self):
         net = op.network.Cubic(shape=[3, 3, 3], connectivity=6)