Repository: pmeal/openpnm
Issue: amalgamate_data fails when key not present

The previous version of amalgamate_data would give the official Python KeyError if the key was not found on network or geometry objects, but is currently return a None Type error or something.  Since @TomTranter broke it he can fix it ;-)


Changed files:
OpenPNM/Utilities/misc.py, OpenPNM/Postprocessing/Plots.py, OpenPNM/Network/__Delaunay__.py, OpenPNM/Base/__Core__.py, OpenPNM/Utilities/vertexops.py, OpenPNM/Algorithms/__OrdinaryPercolation__.py

diff --git a/OpenPNM/Algorithms/__OrdinaryPercolation__.py b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
index 8349663f1..9080ce7c8 100644
--- a/OpenPNM/Algorithms/__OrdinaryPercolation__.py
+++ b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
@@ -53,13 +53,16 @@ class OrdinaryPercolation(GenericAlgorithm):
     To run this algorithm, use 'setup()' to provide the necessary simulation
     """
 
-    def __init__(self,invading_phase=None,defending_phase=None,**kwargs):
+    def __init__(self,invading_phase=None,defending_phase=None,residual_pores=None,residual_throats=None,**kwargs):
         r"""
 
         """
         super(OrdinaryPercolation,self).__init__(**kwargs)
         self._phase_inv = invading_phase
         self._phase_def = defending_phase
+        self._residual_pores = residual_pores
+        self._residual_throats = residual_throats
+        
         logger.debug("Create Drainage Percolation Algorithm Object")
 
     def run(self,
@@ -110,8 +113,10 @@ class OrdinaryPercolation(GenericAlgorithm):
 
         #Create pore and throat conditions lists to store inv_val at which each is invaded
         self._p_inv = sp.zeros((self._net.num_pores(),),dtype=float)
+        self._p_inv.fill(sp.inf)
         self._p_seq = sp.zeros_like(self._p_inv,dtype=int)
         self._t_inv = sp.zeros((self._net.num_throats(),),dtype=float)
+        self._t_inv.fill(sp.inf)
         self._t_seq = sp.zeros_like(self._t_inv,dtype=int)
         #Determine the invasion pressures to apply
         try:
@@ -172,6 +177,8 @@ class OrdinaryPercolation(GenericAlgorithm):
         """
         #Generate a tlist containing boolean values for throat state
         Tinvaded = self._t_cap<=inv_val
+        #if self._residual_throats is not None:
+        #    Tinvaded[self._residual_throats]=True
         #Finding all pores that can be invaded at specified pressure
         clusters = self._net.find_clusters(Tinvaded)
         #Find all pores with at least 1 invaded throat (invaded)
@@ -181,6 +188,8 @@ class OrdinaryPercolation(GenericAlgorithm):
         temp = P12[Tinvaded]
         temp = sp.hstack((temp[:,0],temp[:,1]))
         Pinvaded[temp] = True
+        #if self._residual_pores is not None:
+        #    Pinvaded[self._residual_pores]=True
         if self._AL:
             #Add injection sites to Pinvaded
             Pinvaded[self._inv_sites] = True
@@ -195,12 +204,14 @@ class OrdinaryPercolation(GenericAlgorithm):
             inv_clusters = sp.r_[0:self._net.num_pores()]
         #Store invasion pressure in pores and throats
         pmask = np.in1d(clusters,inv_clusters)
+        #if self._residual_pores is not None:
+        #    pmask[self._residual_pores]==True
         #Store result of invasion step
-        self._p_inv[(self._p_inv==0)*(pmask)] = inv_val
+        self._p_inv[(self._p_inv==sp.inf)*(pmask)] = inv_val
         #Determine Pc_invaded for throats as well
         temp = self._net['throat.conns']
         tmask = (pmask[temp[:,0]] + pmask[temp[:,1]])*(Tinvaded)
-        self._t_inv[(self._t_inv==0)*(tmask)] = inv_val
+        self._t_inv[(self._t_inv==sp.inf)*(tmask)] = inv_val
 
     def evaluate_trapping(self, outlets):
         r"""
@@ -242,7 +253,6 @@ class OrdinaryPercolation(GenericAlgorithm):
                 trapped_throat_array=np.asarray([False]*len(Cstate))
                 trapped_throat_array[trapped_throats]=True
                 self._t_trap[(self._t_trap == 0)*trapped_throat_array] = inv_val
-                #trapped_throats = sp.where(Cstate==2)[0]
                 self._t_trap[(self._t_trap == 0)*(Cstate==2)] = inv_val
         self._p_inv[self._p_trap > 0] = sp.inf
         self._t_inv[self._t_trap > 0] = sp.inf
@@ -336,12 +346,16 @@ class OrdinaryPercolation(GenericAlgorithm):
               Pc = PcPoints[i]
               Snwp_p[i] = sum(Pvol[self._p_inv[pores]<=Pc])/Pvol_tot
               Snwp_t[i] = sum(Tvol[self._t_inv[throats]<=Pc])/Tvol_tot
+          if sp.mean(self._phase_inv["pore.contact_angle"]) < 90:
+              Snwp_p = 1 - Snwp_p
+              Snwp_t = 1 - Snwp_t
+              PcPoints *= -1
           plt.plot(PcPoints,Snwp_p,'r.-')
           plt.plot(PcPoints,Snwp_t,'b.-')
           r'''
           TODO: Add legend to distinguish the pore and throat curves
           '''
-          plt.xlim(xmin=0)
+          #plt.xlim(xmin=0)
           plt.show()
 
 
diff --git a/OpenPNM/Base/__Core__.py b/OpenPNM/Base/__Core__.py
index b06d4eaf9..e11a81f01 100644
--- a/OpenPNM/Base/__Core__.py
+++ b/OpenPNM/Base/__Core__.py
@@ -199,7 +199,7 @@ class Core(Base):
 
         '''
         #Check no duplicate or invalid locations
-        locs = sp.unique(new_order.values())[0]
+        locs = sp.unique(new_order.values())
         if len(locs) < len(new_order.values()):
             raise Exception('Duplicates found in the order')
 
@@ -1042,9 +1042,37 @@ class Core(Base):
         This makes an effort to maintain the data 'type' when possible; however
         when data is missing this can be tricky.  Float and boolean data is
         fine, but missing ints are converted to float when nans are inserted.
+        
+        Examples
+        --------
+        >>> import OpenPNM
+        >>> pn = OpenPNM.Network.TestNet()
+        >>> Ps = pn.pores('top',mode='not')
+        >>> Ts = pn.find_neighbor_throats(pores=Ps,mode='intersection',flatten=True)
+        >>> geom = OpenPNM.Geometry.TestGeometry(network=pn,pores=Ps,throats=Ts)
+        >>> Ps = pn.pores('top')
+        >>> Ts = Ts = pn.find_neighbor_throats(pores=Ps,mode='not_intersection')
+        >>> boun = OpenPNM.Geometry.Boundary(network=pn,pores=Ps,throats=Ts)
+        >>> geom['pore.test_int'] = sp.random.randint(0,100,geom.Np)
+        >>> print(pn['pore.test_int'].dtype)
+        float64
+        >>> boun['pore.test_int'] = sp.ones(boun.Np).astype(int)
+        >>> print(pn['pore.test_int'].dtype)
+        int32
+        >>> boun['pore.test_int'] = sp.rand(boun.Np)<0.5
+        >>> print(pn['pore.test_int'].dtype)
+        bool
+        >>> geom['pore.test_bool'] = sp.rand(geom.Np)<0.5
+        >>> print(pn['pore.test_bool'].dtype)
+        bool
+        >>> boun['pore.test_bool'] = sp.ones(boun.Np).astype(int)
+        >>> print(pn['pore.test_bool'].dtype)
+        bool
+        >>> boun['pore.test_bool'] = sp.rand(boun.Np)<0.5
+        >>> print(pn['pore.test_bool'].dtype)
+        bool
         '''
         element = prop.split('.')[0]
-
         Np = sp.shape(self['pore.all'])[0]
         Nt = sp.shape(self['throat.all'])[0]
         self._temp = {}
@@ -1057,6 +1085,7 @@ class Core(Base):
         dtypes = []
         dtypenames = []
         prop_found = False  #Flag to indicate if prop was found on a sub-object
+        nan_locs = sp.zeros(len(temp)).astype('bool')
         for item in sources:
             #Check if sources were given as list of objects OR names
             try: item.name
@@ -1066,6 +1095,7 @@ class Core(Base):
                 values = sp.ones_like(locations)*sp.nan
                 dtypenames.append('nan')
                 dtypes.append(sp.dtype(bool))
+                nan_locs[locations]=True
             else:
                 values = item[prop]
                 prop_found = True
@@ -1078,23 +1108,36 @@ class Core(Base):
         for i in range(len(temp)):
             cnames.append(temp[i].__class__.__name__)
             shapes.append(sp.shape(temp[i]))
-            if cnames[i] not in ['int','float','bool','int32','int64','float32','float64','NoneType']:
+            if temp[i].__class__.__name__ not in ['int','int32','int64','float','float32','float64','bool','bool_','NoneType']:
                 make_object = True
                 try:
                     dtypes.append(temp[i].dtype.name)
                 except AttributeError:
                     dtypes.append('NoneType')
-        #if all the data was in an acceptable format to return straight away
+        cnames = sp.asarray(cnames)
+        #if all the data was in an acceptable format to return straight away - check data mixtures
         if make_object == False:
+            #All same data type
             if len(sp.unique(cnames))==1:
                 temp_cname = cnames[0]
+            elif len(sp.unique(dtypes))==1:
+                temp_cname = dtypes[0]
+            #If any data is bool - preserve this and convert rest to False
+            elif len(sp.unique(cnames))>1 and sp.logical_or('bool' in cnames,'bool_' in cnames):
+                temp_cname = 'bool'
+                nan_locs[~sp.logical_or(cnames=='bool',cnames=='bool_')]=True
+                logger.warning(prop+' has been converted to bool, some data may be lost')
+            #If mixture with no bool convert to float
             else:
-                #convert a mixture to float
                 temp_cname='float'
+                logger.warning(prop+' has been converted to float.')
         else:
             #keep as object
             temp_cname = 'object'
         temp=temp.astype(temp_cname)
+        if sum(nan_locs) > 0:
+            if sp.logical_or(temp_cname=='bool',temp_cname=='bool_'):
+                temp[nan_locs]=False
         #If we have all the data in ndarray of same type and shape we can re-cast into the proper shape
         if (len(sp.unique(cnames)) == 1)&(cnames[0]=='ndarray')&(len(sp.unique(shapes)) == 1):
             temp_2_shape = list(shapes[0])
diff --git a/OpenPNM/Network/__Delaunay__.py b/OpenPNM/Network/__Delaunay__.py
index bc0f28101..4da4434d1 100644
--- a/OpenPNM/Network/__Delaunay__.py
+++ b/OpenPNM/Network/__Delaunay__.py
@@ -103,12 +103,16 @@ class Delaunay(GenericNetwork):
         #coords = sp.rand(self._Np,3)*[self._Lx,self._Ly,self._Lz]
         #Seeding Code
         coords = np.zeros([self._Np,3])
-        i = 0
-        while i < self._Np:
-            coord = np.array([np.random.uniform(0,self._Lx,1),np.random.uniform(0,self._Ly,1),np.random.uniform(0,self._Lz,1)]).T
-            if self._reject(coord/np.array([self._Lx,self._Ly,self._Lz])) == False:
-                coords[i]=coord
-                i += 1
+        #reject points close to boundaries - if False there will be slightly more
+        rejection = [False,False,True]
+        for j in range(3):
+            i = 0
+            while i < self._Np:
+                coord = np.random.uniform(0,1,1)
+                if self._reject(coord) == rejection[j]:
+                    coords[i][j]=coord
+                    i += 1
+        coords*=np.array([self._Lx,self._Ly,self._Lz])
         #Seeding Code
         #Uniform Random Generator
         #coords = np.array([np.random.uniform(0,self._Lx,self._Np),np.random.uniform(0,self._Ly,self._Np),np.random.uniform(0,self._Lz,self._Np)]).T
@@ -125,15 +129,10 @@ class Delaunay(GenericNetwork):
 
     def _reject(self,point):
 
-        x = point[0,0]
-        y = point[0,1]
-        #z = point[0,2]
-        Px = self._prob_func(x)
-        Py = self._prob_func(y)
-        #Pz = prob_func(z)
+        P = self._prob_func(point)
         nrand = np.random.uniform(0,1,1)
-
-        if Px < nrand and Py < nrand:
+        #place more points at the sides of the domain and fewer at the top and bottom
+        if P < nrand:
             rejection = True
         else:
             rejection = False
@@ -555,13 +554,14 @@ class Delaunay(GenericNetwork):
         --------
         >>> import OpenPNM
         >>> pn = OpenPNM.Network.Delaunay(num_pores=100, domain_size=[3,2,1])
+        >>> import OpenPNM.Utilities.vertexops as vo
         >>> pn.add_boundaries()
         >>> B1 = pn.pores("left_boundary")
         >>> B2 = pn.pores("right_boundary")
         >>> pn.domain_length(B1,B2)
         2.0
         """
-        L = self.vertex_dimension(face_1,face_2,parm='length')
+        L = vo.vertex_dimension(self,face_1,face_2,parm='length')
         return L
 
     def domain_area(self,face):
@@ -578,7 +578,7 @@ class Delaunay(GenericNetwork):
         >>> pn.domain_area(B1)
         3.0
         """
-        A = self.vertex_dimension(face,parm='area')
+        A = vo.vertex_dimension(self,face,parm='area')
 
         return A
 
diff --git a/OpenPNM/Postprocessing/Plots.py b/OpenPNM/Postprocessing/Plots.py
index 5c7637bf7..8ff5ed9b0 100644
--- a/OpenPNM/Postprocessing/Plots.py
+++ b/OpenPNM/Postprocessing/Plots.py
@@ -176,11 +176,11 @@ def drainage_curves(inv_alg,
     >>> phys2 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase2,pores=pn.pores(),throats=pn.throats())
     >>> IP = OpenPNM.Algorithms.InvasionPercolation(network=pn)
     >>> IP.run(invading_phase=phase1, defending_phase=phase2, inlets=pn.pores('top'), outlets=pn.pores('bottom'))
-    IP algorithm at 0 % completion at 0.0 seconds
-    IP algorithm at 20 % completion at 0.0 seconds
-    IP algorithm at 40 % completion at 0.0 seconds
-    IP algorithm at 60 % completion at 0.0 seconds
-    IP algorithm at 100% completion at  0.0  seconds
+         IP algorithm at 0 % completion at 0.0 seconds
+         IP algorithm at 20 % completion at 0.0 seconds
+         IP algorithm at 40 % completion at 0.0 seconds
+         IP algorithm at 60 % completion at 0.0 seconds
+         IP algorithm at 100% completion at  0.0  seconds
     >>> OpenPNM.Postprocessing.Plots.drainage_curves(IP,timing='inv_time')
     """
     inv_throats = inv_alg.toindices(inv_alg['throat.'+seq]>0)
diff --git a/OpenPNM/Utilities/misc.py b/OpenPNM/Utilities/misc.py
index 61ef58a6f..703b847ba 100644
--- a/OpenPNM/Utilities/misc.py
+++ b/OpenPNM/Utilities/misc.py
@@ -186,7 +186,7 @@ def amalgamate_data(objs=[]):
     if type(objs) is not list:
         objs = list(objs)
     data_amalgamated = {}
-    exclusion_list = ['pore.centroid','pore.vertices','throat.centroid','throat.offset_vertices','throat.vertices','throat.normal','throat.perimeter']
+    exclusion_list = ['pore.centroid','pore.vertices','throat.centroid','throat.offset_vertices','throat.vertices','throat.normal','throat.perimeter','pore.vert_index','throat.vert_index']
     for item in objs:
         mro = [module.__name__ for module in item.__class__.__mro__]
         if 'GenericNetwork' in mro: #if Network object, combine Geometry and Network keys
@@ -218,7 +218,7 @@ def amalgamate_data(objs=[]):
                             dict_name = key
                         data_amalgamated.update({dict_name : item[key]})
                 except TypeError:
-                    print(key)
+                    pass
     return data_amalgamated
 
 def conduit_lengths(network,throats=None,mode='pore'):
diff --git a/OpenPNM/Utilities/vertexops.py b/OpenPNM/Utilities/vertexops.py
index 95993b230..bee1f4155 100644
--- a/OpenPNM/Utilities/vertexops.py
+++ b/OpenPNM/Utilities/vertexops.py
@@ -506,13 +506,13 @@ def scale(network,scale_factor=[1,1,1],preserve_vol=True):
     >>> pn.add_boundaries()
     >>> B1 = pn.pores("left_boundary")
     >>> B2 = pn.pores("right_boundary")
-    >>> Vol = pn.vertex_dimension(B1,B2)
+    >>> Vol = vo.vertex_dimension(pn,B1,B2)
     >>> vo.scale(network=pn,scale_factor=[2,1,1])
-    >>> Vol2 = pn.vertex_dimension(B1,B2)
+    >>> Vol2 = vo.vertex_dimension(pn,B1,B2)
     >>> np.around(Vol-Vol2,5)
     0.0
     >>> vo.scale(network=pn,scale_factor=[2,1,1],preserve_vol=False)
-    >>> Vol3 = pn.vertex_dimension(B1,B2)
+    >>> Vol3 = vo.vertex_dimension(pn,B1,B2)
     >>> np.around(Vol3/Vol,5)
     2.0