Repository: PMEAL/OpenPNM
Pull request #649 (2016-09-28 18:49:18)

Fix cubic dual
There was a major bug in this class (the network size was being hard coded, presumably due to some left over debugging lines).  This could have been avoided it we'd done proper team-based code review before merging.  Also, I took the opportunity to adjust how the duals were generated, so now requesting a shape of [4,4,4] returns a [4,4,4].  It was previously being padded to [5,5,5], which was not very logical.


Patch:
diff --git a/OpenPNM/Network/__CubicDual__.py b/OpenPNM/Network/__CubicDual__.py
index 79450f433..81cf3d755 100644
--- a/OpenPNM/Network/__CubicDual__.py
+++ b/OpenPNM/Network/__CubicDual__.py
@@ -6,8 +6,8 @@ CubicDual: Generate a cubic lattice with an interpentrating dual network
 
 """
 import scipy as sp
-import OpenPNM.Utilities.misc as misc
-from OpenPNM.Network import GenericNetwork
+from OpenPNM.Network.tools import stitch
+from OpenPNM.Network import GenericNetwork, Cubic
 from OpenPNM.Base import logging
 logger = logging.getLogger(__name__)
 
@@ -45,51 +45,49 @@ class CubicDual(GenericNetwork):
     def __init__(self, shape=None, spacing=[1, 1, 1], label_1='primary',
                  label_2='secondary', **kwargs):
         super().__init__(**kwargs)
-        import OpenPNM as op
-        spacing = sp.array(1)
-        shape = sp.array([5, 5, 5])
-        label_1 = 'primary'
-        label_2 = 'secondary'
         spacing = sp.array(spacing)
         shape = sp.array(shape)
-        net = op.Network.Cubic(shape=shape, spacing=[1, 1, 1])
-        net.add_boundaries()
+        # Deal with non-3D shape arguments
+        shape = sp.pad(shape, [0, 3-shape.size], mode='constant', constant_values=1)
+        net = Cubic(shape=shape, spacing=[1, 1, 1])
         net['throat.'+label_1] = True
         net['pore.'+label_1] = True
-        dual = op.Network.Cubic(shape=shape+1)
+        single_dim = shape == 1
+        shape[single_dim] = 2
+        dual = Cubic(shape=shape-1, spacing=[1, 1, 1])
+        faces = [['front', 'back'], ['left', 'right'], ['top', 'bottom']]
+        faces = [faces[i] for i in sp.where(~single_dim)[0]]
+        faces = sp.array(faces).flatten().tolist()
+        dual.add_boundaries(faces)
+        # Add secondary network name as a label
         dual['pore.'+label_2] = True
         dual['throat.'+label_2] = True
-        dual['pore.coords'] -= 0.5
-        op.Network.tools.stitch(net, dual, P_network=net.Ps, P_donor=dual.Ps,
-                                len_max=1)
+        # Shift coordinates prior to stitching
+        dual['pore.coords'] += 0.5*(~single_dim)
+        stitch(net, dual, P_network=net.Ps, P_donor=dual.Ps, len_max=1)
         net['throat.interconnect'] = net['throat.stitched']
+        del net['throat.stitched']
         net['pore.coords'] *= spacing
-
-        # Clean-up
+        # Clean-up labels
         net['pore.surface'] = False
         net['throat.surface'] = False
-        surface_labels = ['top', 'bottom', 'front', 'back', 'left', 'right']
-        for face in surface_labels:
-            Ps = net.pores(labels=[face, label_1], mode='intersection')
+        for face in faces:
+            # Remove face label from secondary network since it's internal now
+            Ps = net.pores(labels=[face, label_2], mode='intersection')
             net['pore.'+face][Ps] = False
             Ps = net.pores(labels=[face+'_boundary'])
-            net['pore.surface'][Ps] = True
             net['pore.'+face][Ps] = True
-
-        [net.pop(item) for item in list(net.keys()) if 'boundary' in item]
-        # Label all remaining 'face' pores as 'surface'
-        Ps = net.pores(labels=surface_labels)
-        net['pore.surface'][Ps] = True
-        Ts = net.find_neighbor_throats(pores=net.pores('surface'),
-                                       mode='intersection')
-        net['throat.surface'][Ts] = True
+            Ps = net.pores(face)
+            net['pore.surface'][Ps] = True
+            Ts = net.find_neighbor_throats(pores=Ps, mode='intersection')
+            net['throat.surface'][Ts] = True
+            net['throat.'+face] = net.tomask(throats=Ts)
+        [net.pop(item) for item in net.labels() if 'boundary' in item]
         # Label non-surface pores and throats as internal
         net['pore.internal'] = ~net['pore.surface']
         Ts = net.find_neighbor_throats(pores=net['pore.internal'])
         net['throat.internal'] = False
         net['throat.internal'][Ts] = True
-        # Remove unused labels
-        del net['throat.stitched']
         # Transfer all dictionary items from 'net' to 'self'
         [self.update({item: net[item]}) for item in net]
         del self.workspace[net.name]
diff --git a/OpenPNM/Network/tools.py b/OpenPNM/Network/tools.py
index 04e493185..c0b60e2c7 100644
--- a/OpenPNM/Network/tools.py
+++ b/OpenPNM/Network/tools.py
@@ -345,6 +345,7 @@ def stitch(network, donor, P_network, P_donor, method='nearest',
     # Ensure Networks have no associated objects yet
     if (len(network._simulation()) > 1) or (len(donor._simulation()) > 1):
         raise Exception('Cannot stitch a Network with active sibling objects')
+    network['throat.stitched'] = False
     # Get the initial number of pores and throats
     N_init = {}
     N_init['pore'] = network.Np


Test patch:
diff --git a/test/unit/Network/CubicDualTest.py b/test/unit/Network/CubicDualTest.py
index bdc444551..55419881d 100644
--- a/test/unit/Network/CubicDualTest.py
+++ b/test/unit/Network/CubicDualTest.py
@@ -4,25 +4,83 @@ import scipy as sp
 
 class CubicDualTest:
 
-    def test_generation(self):
+    def test_generation_3D(self):
         net = op.Network.CubicDual(shape=[5, 5, 5], label_1='primary',
                                    label_2='secondary')
-        assert net.Np == 491
-        assert net.Nt == 2590
-        assert net.num_pores('all') == 491
-        assert net.num_pores('back') == 61
-        assert net.num_pores('bottom') == 61
-        assert net.num_pores('front') == 61
-        assert net.num_pores('internal') == 189
-        assert net.num_pores('left') == 61
-        assert net.num_pores('primary') == 275
-        assert net.num_pores('right') == 61
-        assert net.num_pores('secondary') == 216
-        assert net.num_pores('surface') == 302
-        assert net.num_pores('top') == 61
-        assert net.num_throats('all') == 2590
-        assert net.num_throats('interconnect') == 1600
-        assert net.num_throats('internal') == 1690
-        assert net.num_throats('primary') == 450
-        assert net.num_throats('secondary') == 540
-        assert net.num_throats('surface') == 900
+        assert net.Np == 285
+        assert net.Nt == 1436
+        assert net.num_pores('all') == 285
+        assert net.num_pores('back') == 41
+        assert net.num_pores('bottom') == 41
+        assert net.num_pores('front') == 41
+        assert net.num_pores('internal') == 91
+        assert net.num_pores('left') == 41
+        assert net.num_pores('primary') == 125
+        assert net.num_pores('right') == 41
+        assert net.num_pores('secondary') == 160
+        assert net.num_pores('surface') == 194
+        assert net.num_pores('top') == 41
+        assert net.num_throats('all') == 1436
+        assert net.num_throats('interconnect') == 896
+        assert net.num_throats('internal') == 860
+        assert net.num_throats('primary') == 300
+        assert net.num_throats('secondary') == 240
+        assert net.num_throats('surface') == 576
+        assert net.num_throats('top') == 104
+        assert net.num_throats('bottom') == 104
+        assert net.num_throats('left') == 104
+        assert net.num_throats('right') == 104
+        assert net.num_throats('front') == 104
+        assert net.num_throats('back') == 104
+
+    def test_generation_2D_XY(self):
+        net = op.Network.CubicDual(shape=[5, 5, 1], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        assert net.Nt == 176
+        assert net.num_pores('left') == 9
+        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
+        net.num_throats('interconnect') == 96
+
+    def test_generation_2D_XZ(self):
+        net = op.Network.CubicDual(shape=[5, 1, 5], label_1='primary',
+                                   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
+        assert net.num_pores('bottom') == 9
+        net.num_throats('interconnect') == 96
+
+    def test_generation_2D_YZ(self):
+        net = op.Network.CubicDual(shape=[1, 5, 5], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        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
+        net.num_throats('interconnect') == 96
+
+    def test_generation_2D_2_dims(self):
+        net = op.Network.CubicDual(shape=[5, 5], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        assert net.Nt == 176
+        assert net.num_pores('left') == 9
+        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
+        net.num_throats('interconnect') == 96