Repository: PMEAL/OpenPNM
Pull request #1182 (2018-12-27 23:06:49)

fixing performance regression in Cubic


Patch:
diff --git a/openpnm/network/Cubic.py b/openpnm/network/Cubic.py
index 2265abec7..d3e7478a5 100644
--- a/openpnm/network/Cubic.py
+++ b/openpnm/network/Cubic.py
@@ -9,6 +9,8 @@ import numpy as np
 import scipy as sp
 from openpnm.network import GenericNetwork
 from openpnm import topotools
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
 
 
 class Cubic(GenericNetwork):
@@ -93,7 +95,7 @@ class Cubic(GenericNetwork):
         # Store original network shape
         self._shape = sp.shape(arr)
         # Store network spacing
-        spacing = sp.array(spacing)
+        spacing = sp.around(spacing, decimals=15)
         if spacing.size == 2:
             spacing = sp.concatenate((spacing, [1]))
         self._spacing = sp.ones(3)*sp.array(spacing, ndmin=1)
@@ -151,14 +153,26 @@ class Cubic(GenericNetwork):
         self['throat.conns'] = pairs
         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
+        self._label_surface_pores()
+        topotools.label_faces(network=self)
+        Ps = self['pore.surface']
+        self['throat.surface'] = sp.all(Ps[self['throat.conns']], axis=1)
         # Scale network to requested spacing
         self['pore.coords'] *= spacing
 
+    def _label_surface_pores(self):
+        r"""
+        """
+        hits = sp.zeros_like(self.Ps, dtype=bool)
+        dims = topotools.dimensionality(self)
+        mn = sp.amin(self['pore.coords'], axis=0)
+        mx = sp.amax(self['pore.coords'], axis=0)
+        for ax in [0, 1, 2]:
+            if dims[ax]:
+                hits += self['pore.coords'][:, ax] <= mn[ax]
+                hits += self['pore.coords'][:, ax] >= mx[ax]
+        self['pore.surface'] = hits
+
     def add_boundary_pores(self, labels=['top', 'bottom', 'front', 'back',
                                          'left', 'right'], spacing=None):
         r"""
@@ -183,7 +197,7 @@ class Cubic(GenericNetwork):
             labels = [labels]
         x, y, z = self['pore.coords'].T
         if spacing is None:
-            spacing = self._spacing
+            spacing = self._get_spacing()
         else:
             spacing = sp.array(spacing)
             if spacing.size == 1:
@@ -202,36 +216,53 @@ class Cubic(GenericNetwork):
         scale['bottom'] = scale['top'] = [1, 1, 0]
 
         for label in labels:
-            Ps = self.pores(label)
-            topotools.clone_pores(network=self, pores=Ps,
-                                  labels=label+'_boundary')
-            # Translate cloned pores
-            ind = self.pores(label+'_boundary')
-            coords = self['pore.coords'][ind]
-            coords = coords*scale[label] + offset[label]
-            self['pore.coords'][ind] = coords
+            try:
+                Ps = self.pores(label)
+                topotools.clone_pores(network=self, pores=Ps,
+                                      labels=label+'_boundary')
+                # Translate cloned pores
+                ind = self.pores(label+'_boundary')
+                coords = self['pore.coords'][ind]
+                coords = coords*scale[label] + offset[label]
+                self['pore.coords'][ind] = coords
+            except KeyError:
+                logger.warning('No pores labelled ' + label + ' were found, ' +
+                               'skipping boundary addition')
 
     def _get_spacing(self):
         # Find Network spacing
         P12 = self['throat.conns']
         C12 = self['pore.coords'][P12]
-        V = np.abs(np.squeeze(np.diff(C12, axis=1)))
-        if np.any(np.sum(V == 0, axis=1) != 2):
-            raise Exception('A unique value of spacing could not be found')
-        spacing = [None, None, None]
-        for axis in [0, 1, 2]:
-            temp = np.unique(np.around(V[:, axis], decimals=10))
-            if np.size(temp) > 2:
-                raise Exception('A unique value of spacing could not be found')
-            else:
-                spacing[axis] = temp[1]
+        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)
+        spacing = [0, 0, 0]
+        dims = topotools.dimensionality(self)
+        # Ensure vectors point in n-dims unique directions
+        c = {tuple(row): 1 for row in 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)
+        for ax in [0, 1, 2]:
+            if dims[ax]:
+                inds = sp.where(vec[:, ax] == vec[:, ax].max())[0]
+                temp = sp.unique(mag[inds])
+                if np.size(temp) > 1:
+                    raise Exception('A unique value of spacing could not be found')
+                else:
+                    spacing[ax] = temp[0]
         return sp.array(spacing)
 
     spacing = property(fget=_get_spacing)
 
     def _get_shape(self):
         L = np.ptp(self['pore.coords'], axis=0)
-        return L/self.spacing + 1
+        mask = L.astype(bool)
+        S = self.spacing
+        shape = sp.array([1, 1, 1], int)
+        shape[mask] = L[mask]/S[mask] + 1
+        return shape
 
     shape = property(fget=_get_shape)
 
diff --git a/openpnm/topotools/__init__.py b/openpnm/topotools/__init__.py
index d81cfe440..60d1b44bb 100644
--- a/openpnm/topotools/__init__.py
+++ b/openpnm/topotools/__init__.py
@@ -42,6 +42,8 @@ from .topotools import plot_networkx
 from .topotools import reduce_coordination
 from .topotools import reflect_base_points
 from .topotools import remove_isolated_clusters
+from .topotools import rotate_coords
+from .topotools import shear_coords
 from .topotools import site_percolation
 from .topotools import stitch
 from .topotools import subdivide
diff --git a/openpnm/topotools/topotools.py b/openpnm/topotools/topotools.py
index 5453b3040..4e9094b4c 100644
--- a/openpnm/topotools/topotools.py
+++ b/openpnm/topotools/topotools.py
@@ -807,6 +807,117 @@ def bond_percolation(ij, occupied_bonds):
     return tup(s_labels, b_labels)
 
 
+def rotate_coords(network, a=0, b=0, c=0, R=None):
+    r"""
+    Rotates coordinates a given amount about each axis
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network whose pore coordinates should be transformed
+    a : scalar
+        The amount in degrees to rotate about the x-axis
+    b : scalar
+        The amount in degrees to rotate about the y-axis
+    c : scalar
+        The amount in degrees to rotate about the z-axis
+    R : array_like
+        Rotation matrix.  Must be a 3-by-3 matrix since pore coordinates are
+        always in 3D.  If this is given then the other individual arguments
+        are ignored.
+
+    See Also
+    --------
+    rotate_coords
+
+    Notes
+    -----
+    It is possible to rotate about anyof the three axes by specifying ``a``,
+    ``b``, and/or ``c``.  In this case each rotation is applied in sequence.
+
+    """
+    if R is None:
+        if a:
+            R = sp.array([[1, 0, 0],
+                          [0, sp.cos(sp.deg2rad(a)), -sp.sin(sp.deg2rad(a))],
+                          [0, sp.sin(sp.deg2rad(a)), sp.cos(sp.deg2rad(a))]])
+            network['pore.coords'] = sp.tensordot(network['pore.coords'], R,
+                                                  axes=(1, 1))
+        if b:
+            R = sp.array([[sp.cos(sp.deg2rad(b)), 0, -sp.sin(sp.deg2rad(b))],
+                          [0, 1, 0],
+                          [sp.sin(sp.deg2rad(b)), 0, sp.cos(sp.deg2rad(b))]])
+            network['pore.coords'] = sp.tensordot(network['pore.coords'], R,
+                                                  axes=(1, 1))
+        if c:
+            R = sp.array([[sp.cos(sp.deg2rad(c)), -sp.sin(sp.deg2rad(c)), 0],
+                          [sp.sin(sp.deg2rad(c)), sp.cos(sp.deg2rad(c)), 0],
+                          [0, 0, 1]])
+            network['pore.coords'] = sp.tensordot(network['pore.coords'], R,
+                                                  axes=(1, 1))
+    else:
+        network['pore.coords'] = sp.tensordot(network['pore.coords'], R,
+                                              axes=(1, 1))
+
+
+def shear_coords(network, ay=0, az=0, bx=0, bz=0, cx=0, cy=0, S=None):
+    r"""
+    Shears the coordinates a given amount about along axis
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network whose pore coordinates should be transformed
+    ay : scalar
+        The factor by which to shear along the x-axis as a function of y
+    az : scalar
+        The factor by which to shear along the x-axis as a function of z
+    bx : scalar
+        The factor by which to shear along the y-axis as a function of x
+    bz : scalar
+        The factor by which to shear along the y-axis as a function of z
+    cx : scalar
+        The factor by which to shear along the z-axis  as a function of x
+    cy : scalar
+        The factor by which to shear along the z-axis as a function of y
+    S : array_like
+        The shear matrix.  Must be a 3-by-3 matrix since pore coordinates are
+        always in 3D.  If this is given then the other individual arguments
+        are ingnored.
+
+    See Also
+    --------
+    rotate_coords
+
+    Notes
+    -----
+    The shear along the i *th* -axis is given as i\* = i + aj.  This means
+    the new i coordinate is the old one plus some linear factor *a* in the
+    j *th* direction.
+
+    The values of ``a``, ``b``, and ``c`` are essentially the inverse of the
+    slope to be formed by the neighboring layers of sheared pores.  A value of
+    0 means no shear, and neighboring points are stacked directly on top of
+    each other; a value of 1 means they form a 45 degree diagonal, and so on.
+
+    If ``S`` is given, then is should be of the form:
+
+    ::
+
+        S = [[1 , ay, az],
+             [bx, 1 , bz],
+             [cx, cy, 1 ]]
+
+        where any of the off-diagonal components can be 0
+
+    """
+    if S is None:
+        S = sp.array([[1, ay, az],
+                      [bx, 1, bz],
+                      [cx, cy, 1]])
+    network['pore.coords'] = (S@network['pore.coords'].T).T
+
+
 def trim(network, pores=[], throats=[]):
     '''
     Remove pores or throats from the network.


Test patch:
diff --git a/tests/unit/models/geometry/ThroatEndpointsTest.py b/tests/unit/models/geometry/ThroatEndpointsTest.py
index b66dd6b3d..178b2e184 100644
--- a/tests/unit/models/geometry/ThroatEndpointsTest.py
+++ b/tests/unit/models/geometry/ThroatEndpointsTest.py
@@ -350,7 +350,7 @@ class ThroatEndpointsTest:
                                               domain_length=length)
             conds.append(D[0])
             nets.append(net)
-        assert np.all(np.around(np.asarray(conds), 12) == D_ab)
+        assert np.allclose(np.asarray(conds), D_ab)
 
 
 if __name__ == '__main__':
diff --git a/tests/unit/network/CubicTest.py b/tests/unit/network/CubicTest.py
new file mode 100644
index 000000000..aff91fff8
--- /dev/null
+++ b/tests/unit/network/CubicTest.py
@@ -0,0 +1,124 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class CubicTest:
+    def setup_class(self):
+        pass
+
+    def teardown_class(self):
+        pass
+
+    def test_spacing_2D(self):
+        net = op.network.Cubic(shape=[5, 5, 1], spacing=[1, 1])
+        assert sp.all(net.spacing == [1.0, 1.0, 0.0])
+
+    def test_spacing_3D(self):
+        net = op.network.Cubic(shape=[5, 5, 5], spacing=[1, 1, 1])
+        assert sp.all(net.spacing == [1.0, 1.0, 1.0])
+
+    def test_spacing_2D_uneven(self):
+        net = op.network.Cubic(shape=[5, 5, 1], spacing=[1, 2])
+        assert sp.all(net.spacing == [1.0, 2.0, 0.0])
+
+    def test_spacing_3D_uneven(self):
+        net = op.network.Cubic(shape=[3, 4, 5], spacing=[1, 2, 3])
+        assert sp.all(net.spacing == [1.0, 2.0, 3.0])
+
+    def test_shape_1D(self):
+        net = op.network.Cubic(shape=[5, 1, 1])
+        assert sp.all(net.shape == [5, 1, 1])
+        net = op.network.Cubic(shape=[1, 1, 5])
+        assert sp.all(net.shape == [1, 1, 5])
+        net = op.network.Cubic(shape=[1, 5, 1])
+        assert sp.all(net.shape == [1, 5, 1])
+
+    def test_shape_2D(self):
+        net = op.network.Cubic(shape=[5, 5, 1])
+        assert sp.all(net.shape == [5, 5, 1])
+        net = op.network.Cubic(shape=[5, 1, 5])
+        assert sp.all(net.shape == [5, 1, 5])
+        net = op.network.Cubic(shape=[1, 5, 5])
+        assert sp.all(net.shape == [1, 5, 5])
+
+    def test_shape_3D(self):
+        net = op.network.Cubic(shape=[5, 5, 5])
+        assert sp.all(net.shape == [5, 5, 5])
+        net = op.network.Cubic(shape=[3, 4, 5])
+        assert sp.all(net.shape == [3, 4, 5])
+        net = op.network.Cubic(shape=[1, 5, 1])
+        assert sp.all(net.shape == [1, 5, 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])
+
+    def test_spacing_3D_rotated_uneven(self):
+        net = op.network.Cubic(shape=[3, 4, 5], spacing=[1, 2, 3])
+        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, 2.0, 3.0])
+
+    def test_spacing_2D_sheared(self):
+        net = op.network.Cubic(shape=[5, 5, 1], spacing=1)
+        S = sp.array([[1, 1, 0],
+                      [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])
+
+    def test_spacing_2D_sheared_uneven(self):
+        net = op.network.Cubic(shape=[5, 5, 1], spacing=[1, 2])
+        S = sp.array([[1, 1, 0],
+                      [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), 0.0])
+
+    def test_spacing_3D_sheared(self):
+        net = op.network.Cubic(shape=[5, 5, 3], spacing=1)
+        S = sp.array([[1, 1, 0],
+                      [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, 1.0])
+
+    def test_spacing_3D_sheared_uneven(self):
+        net = op.network.Cubic(shape=[3, 4, 5], spacing=[1, 2, 3])
+        S = sp.array([[1, 1, 0],
+                      [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])
+
+    def test_spacing_on_joggled_network(self):
+        net = op.network.Cubic(shape=[3, 4, 5])
+        net['pore.coords'] += sp.rand(net.Np, 3)
+        with pytest.raises(Exception):
+            net.spacing
+
+    def test_spacing_on_network_with_boundary_pores(self):
+        net = op.network.Cubic(shape=[3, 4, 5])
+        net.add_boundary_pores()
+        with pytest.raises(Exception):
+            net.spacing
+
+
+if __name__ == '__main__':
+
+    t = CubicTest()
+    t.setup_class()
+    self = t
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/topotools/TransformationsTest.py b/tests/unit/topotools/TransformationsTest.py
new file mode 100644
index 000000000..6b941fcc8
--- /dev/null
+++ b/tests/unit/topotools/TransformationsTest.py
@@ -0,0 +1,59 @@
+import openpnm as op
+import numpy as np
+from numpy.testing import assert_allclose
+from openpnm import topotools
+
+
+class TransformationsTest:
+
+    def setup_class(self):
+        self.ws = op.Workspace()
+
+    def teardown_class(self):
+        self.ws.clear()
+
+    def test_rotate_points(self):
+        net = op.network.Cubic(shape=[3, 3, 3])
+        topotools.rotate_coords(network=net, a=45)
+        a = np.array([[5.00000000e-01,  0.00000000e+00, 7.07106781e-01],
+                      [5.00000000e-01, -7.07106781e-01, 1.41421356e+00]])
+        assert np.allclose(net['pore.coords'][0:2, :], a)
+        net = op.network.Cubic(shape=[3, 3, 3])
+        topotools.rotate_coords(network=net, a=45, b=45, c=45)
+        b = np.array([[-0.10355339, -0.10355339, 0.85355339],
+                      [+0.04289322, -0.95710678, 1.35355339]])
+        assert np.allclose(net['pore.coords'][0:2, :], b)
+        net = op.network.Cubic(shape=[3, 3, 3])
+        q = 45
+        R = np.array([[1, 0, 0],
+                      [0, np.cos(np.deg2rad(q)), -np.sin(np.deg2rad(q))],
+                      [0, np.sin(np.deg2rad(q)), np.cos(np.deg2rad(q))]])
+        topotools.rotate_coords(network=net, R=R)
+        assert np.allclose(net['pore.coords'][0:2, :], a)
+
+    def test_shear_points(self):
+        net = op.network.Cubic(shape=[3, 3, 3])
+        a = np.array([[0.5, 0.5, 0.5],
+                      [0.5, 0.5, 1.5]])
+        assert np.allclose(net['pore.coords'][0:2, :], a)
+        topotools.shear_coords(network=net, ay=1)
+        b = np.array([[1.0, 0.5, 0.5],
+                      [1.0, 0.5, 1.5]])
+        assert np.allclose(net['pore.coords'][0:2, :], b)
+        net = op.network.Cubic(shape=[3, 3, 3])
+        S = [[1.0, 1.0, 0.0],
+             [0.0, 1.0, 0.0],
+             [0.0, 0.0, 1.0]]
+        topotools.shear_coords(network=net, S=S)
+        assert np.allclose(net['pore.coords'][0:2, :], b)
+
+
+if __name__ == '__main__':
+
+    t = TransformationsTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()