Repository: pmeal/openpnm
Issue: Interleave data could infer the wrong shape

In line 803, depending on which object comes first, the inferred shape could be wrong. For instance, if we have `(Nt,1)` and `(Nt,2)` conductances, it'll break (the opposite accidentally works fine, since it'll broadcast the values). See #1676.

https://github.com/PMEAL/OpenPNM/blob/8fbac00f5a01e6b4b57f74d6c89d00a64a594a3d/openpnm/core/Base.py#L800-L805

Changed files:
tests/unit/core/InterleaveDataTest.py, openpnm/core/_base.py

diff --git a/openpnm/core/_base.py b/openpnm/core/_base.py
index 3b64fb4a5..1a4a429b3 100644
--- a/openpnm/core/_base.py
+++ b/openpnm/core/_base.py
@@ -682,14 +682,14 @@ class Base(dict):
         [False False False False]
 
         """
-        # Fetch sources list depending on type of self
+        # Fetch subdomains list depending on type of self
         proj = self.project
         if self._isa() in ['network', 'geometry']:
-            sources = list(proj.geometries().values())
+            subdomains = list(proj.geometries().values())
         elif self._isa() in ['phase', 'physics']:
-            sources = list(proj.find_physics(phase=self))
+            subdomains= list(proj.find_physics(phase=self))
         elif self._isa() in ['algorithm', 'base']:
-            sources = [self]
+            subdomains= [self]
         else:
             raise Exception('Unrecognized object type, cannot find dependents')
 
@@ -698,23 +698,26 @@ class Base(dict):
         N = self.project.network._count(element)
 
         # Attempt to fetch the requested array from each object
-        arrs = [obj.get(prop, None) for obj in sources]
+        arrs = [obj.get(prop, None) for obj in subdomains]
 
         # Check for missing sources, and add None to arrs if necessary
-        if N > sum([obj._count(element) for obj in sources]):
+        if N > sum([obj._count(element) for obj in subdomains]):
             arrs.append(None)
 
         # Obtain list of locations for inserting values
-        locs = [self._get_indices(element, item.name) for item in sources]
+        locs = [self._get_indices(element, item.name) for item in subdomains]
 
         if np.all([item is None for item in arrs]):  # prop not found anywhere
             raise KeyError(prop)
 
         # Let's start by handling the easy cases first
         if not any([a is None for a in arrs]):
-            # All objs present and array found on all objs
-            shape = list(arrs[0].shape)
-            shape[0] = N
+            # All objs are present and array found on all objs
+            try:
+                W = max([a.shape[1] for a in arrs])  # Width of array
+                shape = [N, W]
+            except:
+                shape = [N, ]
             types = [a.dtype for a in arrs]
             if len(set(types)) == 1:
                 # All types are the same
diff --git a/tests/unit/core/InterleaveDataTest.py b/tests/unit/core/InterleaveDataTest.py
index ad0d57d1e..83f521919 100644
--- a/tests/unit/core/InterleaveDataTest.py
+++ b/tests/unit/core/InterleaveDataTest.py
@@ -301,6 +301,29 @@ class InterleaveDataTest:
         with pytest.raises(Exception):
             pn['pore.bee.bop'] = 1
 
+    def test_interleaving_nt_by_1_with_nt_by_2(self):
+        pn = op.network.Cubic(shape=[10, 1, 1])
+        geo1 = op.geometry.SpheresAndCylinders(network=pn, pores=range(5),
+                                               throats=range(4))
+        geo2 = op.geometry.SpheresAndCylinders(network=pn, pores=range(5, 10),
+                                               throats=range(4, 9))
+        geo1['throat.conduit_data'] = np.ones((geo1.Nt, 2))
+        geo2['throat.conduit_data'] = np.ones((geo2.Nt, 1))*2
+        a = np.array([[1., 1.],
+                      [1., 1.],
+                      [1., 1.],
+                      [1., 1.],
+                      [2., 2.],
+                      [2., 2.],
+                      [2., 2.],
+                      [2., 2.],
+                      [2., 2.]])
+        assert np.all(pn['throat.conduit_data'] == a)
+        # Now do it the other way
+        geo1['throat.conduit_data'] = np.ones((geo1.Nt, 1))
+        geo2['throat.conduit_data'] = np.ones((geo2.Nt, 2))*2
+        assert np.all(pn['throat.conduit_data'] == a)
+
 
 if __name__ == '__main__':