Repository: PMEAL/OpenPNM
Pull request #1258 (2019-04-25 23:42:19)

Fix bug in _get_shape, which happens with small spacing values (<1e-9)
Addresses issues #1256 ,#1257.


Patch:
diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
index 37713d66e..90bce2e20 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -746,8 +746,8 @@ class GenericTransport(GenericAlgorithm):
         if not iscoplanar(outlets):
             logger.error('Detected inlet pores are not coplanar')
         tree = cKDTree(data=inlets)
-        Ls = np.unique(np.around(tree.query(x=outlets)[0], decimals=5))
-        if np.size(Ls) != 1:
+        Ls = np.unique(np.float64(tree.query(x=outlets)[0]))
+        if not np.allclose(Ls, Ls[0]):
             logger.error('A unique value of length could not be found')
         length = Ls[0]
         return length
diff --git a/openpnm/network/Cubic.py b/openpnm/network/Cubic.py
index 344fadd30..30ed5ed33 100644
--- a/openpnm/network/Cubic.py
+++ b/openpnm/network/Cubic.py
@@ -90,15 +90,19 @@ class Cubic(GenericNetwork):
     def __init__(self, shape, spacing=[1, 1, 1], connectivity=6, name=None,
                  project=None):
 
+        # Take care of 1D/2D networks
+        shape = sp.array(shape, ndmin=1)
+        shape = np.concatenate((shape, [1]*(3-shape.size))).astype(int)
+
         arr = np.atleast_3d(np.empty(shape))
 
         # Store original network shape
         self._shape = sp.shape(arr)
         # Store network spacing
-        spacing = sp.around(spacing, decimals=15)
+        spacing = sp.float64(spacing)
         if spacing.size == 2:
             spacing = sp.concatenate((spacing, [1]))
-        self._spacing = sp.ones(3)*sp.array(spacing, ndmin=1)
+        self._spacing = sp.ones(3, dtype=float)*sp.array(spacing, ndmin=1)
 
         z = np.tile(np.arange(shape[2]), shape[0]*shape[1])
         y = np.tile(np.repeat(np.arange(shape[1]), shape[2]), shape[0])
@@ -234,22 +238,21 @@ class Cubic(GenericNetwork):
         # Find Network spacing
         P12 = self['throat.conns']
         C12 = self['pore.coords'][P12]
-        mag = np.sqrt(np.sum(np.diff(C12, axis=1)**2, axis=2))
-        vec = sp.around(sp.squeeze(np.diff(C12, axis=1))/mag,
-                        decimals=10)
+        mag = np.linalg.norm(np.diff(C12, axis=1), axis=2)
+        unit_vec = sp.around(sp.squeeze(np.diff(C12, axis=1))/mag, decimals=14)
         spacing = [0, 0, 0]
         dims = topotools.dimensionality(self)
         # Ensure vectors point in n-dims unique directions
-        c = {tuple(row): 1 for row in vec}
+        c = {tuple(row): 1 for row in unit_vec}
         if len(c.keys()) > sum(dims):
             raise Exception('Spacing is undefined when throats point in ' +
                             'more directions than network has dimensions')
-        mag = sp.around(mag.squeeze(), decimals=10)
+        mag = sp.float64(mag.squeeze())
         for ax in [0, 1, 2]:
             if dims[ax]:
-                inds = sp.where(vec[:, ax] == vec[:, ax].max())[0]
+                inds = sp.where(unit_vec[:, ax] == unit_vec[:, ax].max())[0]
                 temp = sp.unique(mag[inds])
-                if np.size(temp) > 1:
+                if not sp.allclose(temp, temp[0]):
                     raise Exception('A unique value of spacing could not be found')
                 else:
                     spacing[ax] = temp[0]
diff --git a/openpnm/topotools/topotools.py b/openpnm/topotools/topotools.py
index 65072ea92..166d8478a 100644
--- a/openpnm/topotools/topotools.py
+++ b/openpnm/topotools/topotools.py
@@ -1311,11 +1311,9 @@ def dimensionality(network):
     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
+    xyz = network["pore.coords"]
+    xyz_unique = [sp.unique(xyz[:, i]) for i in range(3)]
+    return [elem.size != 1 for elem in xyz_unique]
 
 
 def clone_pores(network, pores, labels=['clone'], mode='parents'):


Test patch:
diff --git a/tests/unit/network/CubicTest.py b/tests/unit/network/CubicTest.py
index aff91fff8..8b265fb5f 100644
--- a/tests/unit/network/CubicTest.py
+++ b/tests/unit/network/CubicTest.py
@@ -50,14 +50,18 @@ class CubicTest:
         net = op.network.Cubic(shape=[1, 5, 1])
         assert sp.all(net.shape == [1, 5, 1])
 
+    def test_shape_extremely_small_spacing(self):
+        net = op.network.Cubic(shape=[200, 200, 1], spacing=1e-9)
+        assert sp.allclose(net.shape, [200, 200, 1])
+
     def test_spacing_3D_rotated(self):
         net = op.network.Cubic(shape=[5, 5, 5], spacing=[1, 1, 1])
         theta = 0.1
         R = sp.array([[1, 0, 0],
                       [0, sp.cos(theta), -sp.sin(theta)],
-                      [0, sp.sin(theta), sp.cos(theta)]])
-        net['pore.coords'] = sp.tensordot(net['pore.coords'], R, axes=(1, 1))
-        assert sp.all(net.spacing == [1.0, 1.0, 1.0])
+                      [0, sp.sin(theta), sp.cos(theta)]], dtype=float)
+        net['pore.coords'] = net['pore.coords'] @ R.T
+        assert sp.allclose(net.spacing, [1, 1, 1])
 
     def test_spacing_3D_rotated_uneven(self):
         net = op.network.Cubic(shape=[3, 4, 5], spacing=[1, 2, 3])
@@ -66,7 +70,7 @@ class CubicTest:
                       [0, sp.cos(theta), -sp.sin(theta)],
                       [0, sp.sin(theta), sp.cos(theta)]])
         net['pore.coords'] = sp.tensordot(net['pore.coords'], R, axes=(1, 1))
-        assert sp.all(net.spacing == [1.0, 2.0, 3.0])
+        assert sp.allclose(net.spacing, [1, 2, 3])
 
     def test_spacing_2D_sheared(self):
         net = op.network.Cubic(shape=[5, 5, 1], spacing=1)
@@ -74,7 +78,7 @@ class CubicTest:
                       [0, 1, 0],
                       [0, 0, 1]])
         net['pore.coords'] = (S@net['pore.coords'].T).T
-        assert sp.allclose(net.spacing, [1.0, 2**0.5, 0.0])
+        assert sp.allclose(net.spacing, [1, 2**0.5, 0])
 
     def test_spacing_2D_sheared_uneven(self):
         net = op.network.Cubic(shape=[5, 5, 1], spacing=[1, 2])
@@ -98,7 +102,7 @@ class CubicTest:
                       [0, 1, 0],
                       [0, 0, 1]])
         net['pore.coords'] = (S@net['pore.coords'].T).T
-        assert sp.allclose(net.spacing, [1.0, 2*(2**0.5), 3.0])
+        assert sp.allclose(net.spacing, [1, 2*(2**0.5), 3])
 
     def test_spacing_on_joggled_network(self):
         net = op.network.Cubic(shape=[3, 4, 5])