Repository: PMEAL/OpenPNM
Pull request #1130 (2018-10-16 15:54:15)

Further improvements to merge_pores and connect_pores
`merge_pores` and `connect_pores` now support batch mode, i.e., you can pass a list of list of pores. The reason for doing this is because if you want to consecutively merge/connect a bunch of pores, the adjacency matrix is being rebuilt between each call to `merge_pores` and `connect_pores`. The old API still works, and here's how the new API works:

## Merge pores in batch mode
```
import openpnm as op
network = op.network.Cubic(shape[5, 5, 5])
pores = [[0, 1, 2, 5], [5, 10, 2, 11, 23]]
op.topotools.merge_pores(network=net, pores=pores)
```
After running the above snippet, both the batch of pores `[0, 1, 2, 5]` and `[5, 10, 2, 11, 23]` are merged into two new pores. Note that these two batches share some pores, and ideally the two new pores should be connected because of that. The improved `merge_pores` accounts for this.

## Connect pores in batch mode
```
import openpnm as op
network = op.network.Cubic(shape[5, 5, 5])
ps1 = [[0, 1], [2, 3]]
ps2 = [[5], [7, 9]]
op.topotools.connect_pores(pores1=ps1, pores2=ps2)
```
Running the above snippet results in the following new connections:
```
[0 , 5]     [2 , 7]    [3 , 7]
[1 , 5]     [2 , 9]    [3 , 9]
```

Finally, for the batch mode to work properly, the input pores must be a list of lists or a list of ndarrays, etc. The method throws an exception if input is not properly formatted.


Patch:
diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
index ff996d986..fa4680193 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -524,10 +524,8 @@ class GenericTransport(GenericAlgorithm):
                 import pyamg
             else:
                 raise Exception('pyamg is not installed.')
-            B = np.ones((self.A.shape[0], 1))
-            mc = min(10, int(A.shape[0]//4))
-            ml = pyamg.smoothed_aggregation_solver(A, B, max_coarse=mc)
-            x = ml.solve(b=b, tol=atol)
+            ml = pyamg.ruge_stuben_solver(A)
+            x = ml.solve(b=b, tol=1e-6)
             return x
 
     def results(self, times='all', t_precision=12, **kwargs):
diff --git a/openpnm/topotools/topotools.py b/openpnm/topotools/topotools.py
index 6384e229a..798775b2f 100644
--- a/openpnm/topotools/topotools.py
+++ b/openpnm/topotools/topotools.py
@@ -854,12 +854,15 @@ def trim(network, pores=[], throats=[]):
 #        temp = sp.in1d(network['throat.conns'].flatten(), pores)
 #        temp = sp.reshape(temp, (network.Nt, 2))
 #        Ts = sp.any(temp, axis=1)
+#        Ts = network.Ts[Ts]
+#        tic()
         Ts = network.find_neighbor_throats(pores=~Pkeep, mode='union')
+#        toc()
         if len(Ts) > 0:
             Tkeep[Ts] = False
     if sp.size(throats) > 0:
         Tkeep[throats] = False
-        # The folling IF catches the special case of deleting ALL throats
+        # The following IF catches the special case of deleting ALL throats
         # It removes all throat props, adds 'all', and skips rest of function
         if not sp.any(Tkeep):
             logger.info('Removing ALL throats from network')
@@ -1382,6 +1385,8 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
     Returns the possible connections between two group of pores, and optionally
     makes the connections.
 
+    See ``Notes`` for advanced usage.
+
     Parameters
     ----------
     network : OpenPNM Network Object
@@ -1403,7 +1408,17 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
 
     Notes
     -----
-    It creates the connections in a format which is acceptable by
+    (1) The method also works if ``pores1`` and ``pores2`` are list of lists,
+    in which case it consecutively connects corresponding members of the two
+    lists in a 1-to-1 fashion. Example: pores1 = [[0, 1], [2, 3]] and
+    pores2 = [[5], [7, 9]] leads to creation of the following connections:
+        0 --> 5     2 --> 7     3 --> 7
+        1 --> 5     2 --> 9     3 --> 9
+
+    (2) If you want to use the batch functionality, make sure that each element
+    within ``pores1`` and ``pores2`` are of type list or ndarray.
+
+    (3) It creates the connections in a format which is acceptable by
     the default OpenPNM connection ('throat.conns') and either adds them to
     the network or returns them.
 
@@ -1426,11 +1441,27 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
            [32, 68]])
 
     '''
-    size1 = sp.size(pores1)
-    size2 = sp.size(pores2)
-    array1 = sp.repeat(pores1, size2)
-    array2 = sp.tile(pores2, size1)
-    conns = sp.vstack([array1, array2]).T
+    # Assert that `pores1` and `pores2` are list of lists
+    try:
+        len(pores1[0])
+    except (TypeError, IndexError):
+        pores1 = [pores1]
+    try:
+        len(pores2[0])
+    except (TypeError, IndexError):
+        pores2 = [pores2]
+
+    if len(pores1) != len(pores2):
+        raise Exception('Running in batch mode! pores1 and pores2 must be' + \
+                        ' of the same length.')
+
+    arr1, arr2 = [], []
+    for ps1, ps2 in zip(pores1, pores2):
+        size1 = sp.size(ps1)
+        size2 = sp.size(ps2)
+        arr1.append(sp.repeat(ps1, size2))
+        arr2.append(sp.tile(ps2, size1))
+    conns = sp.vstack([sp.concatenate(arr1), sp.concatenate(arr2)]).T
     if add_conns:
         extend(network=network, throat_conns=conns, labels=labels)
     else:
@@ -1676,20 +1707,40 @@ def merge_pores(network, pores, labels=['merged']):
     # Assert that `pores` is list of lists
     try:
         len(pores[0])
-    except TypeError:
+    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::]
-    for P, NB in zip(Pnew, NBs):
-        connect_pores(network, pores1=P, pores2=NB, labels=labels)
+    
+    # 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 
+    # 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 = [], []    
+    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
+    connect_pores(network, pores2=sp.split(Pnew, N), pores1=NBs, labels=labels)
+    # Trim merged pores from the network
     trim(network=network, pores=sp.concatenate(pores))
 
 
@@ -1847,6 +1898,7 @@ def plot_connections(network, throats=None, fig=None, **kwargs):
 
     """
     import matplotlib.pyplot as plt
+    from mpl_toolkits.mplot3d import Axes3D
 
     if throats is None:
         Ts = network.Ts
@@ -1944,6 +1996,7 @@ def plot_coordinates(network, pores=None, fig=None, **kwargs):
 
     """
     import matplotlib.pyplot as plt
+    from mpl_toolkits.mplot3d import Axes3D
 
     if pores is None:
         Ps = network.Ps


Test patch:
diff --git a/tests/unit/topotools/TopotoolsTest.py b/tests/unit/topotools/TopotoolsTest.py
index bb1e4e67e..067423ca3 100644
--- a/tests/unit/topotools/TopotoolsTest.py
+++ b/tests/unit/topotools/TopotoolsTest.py
@@ -134,6 +134,24 @@ class TopotoolsTest:
         assert_allclose(xyz1, xyz1_desired)
         assert_allclose(xyz2, xyz2_desired)
 
+    def test_connect_pores(self):
+        testnet = op.network.Cubic(shape=[10, 10, 10])
+        Nt_old= testnet.Nt
+        ps1 = [[0, 1], [23, 65]]
+        ps2 = [[55], [982, 555]]
+        topotools.connect_pores(testnet, pores1=ps1, pores2=ps2)
+        am = testnet.create_adjacency_matrix(weights=np.ones(testnet.Nt,
+                                                             dtype=int),
+                                             fmt='csr')
+        conns = testnet['throat.conns']
+        assert len(conns) == Nt_old + 6
+        assert am[0, 55] == 1
+        assert am[1, 55] == 1
+        assert am[23, 982] == 1
+        assert am[23, 555] == 1
+        assert am[65, 982] == 1
+        assert am[65, 555] == 1
+
     def test_ispercolating(self):
         net = op.network.Cubic(shape=[10, 10, 10], connectivity=26)
         tmask = net['throat.all']