Repository: PMEAL/OpenPNM
Pull request #393 (2015-06-16 14:48:28)

Major changes to Voronoi with new voxel based method for calculating pore volumes


Patch:
diff --git a/OpenPNM/Algorithms/__InvasionPercolationTimed__.py b/OpenPNM/Algorithms/__InvasionPercolationTimed__.py
index 41574a33e..82e30ae61 100644
--- a/OpenPNM/Algorithms/__InvasionPercolationTimed__.py
+++ b/OpenPNM/Algorithms/__InvasionPercolationTimed__.py
@@ -306,14 +306,14 @@ class InvasionPercolationTimed(GenericAlgorithm):
         self._plast = len(np.nonzero(self['pore.cluster_final'])[0])
         if self._timing:
             # determine the cluster with the earliest Haines time
-            self._current_cluster = 1 + self._cluster_data['haines_time'].tolist().index(min(self._cluster_data['haines_time']))
+            self._current_cluster = 1 + self._cluster_data['haines_time'].tolist().index(np.min(self._cluster_data['haines_time']))
             # update simulation clock
             logger.debug( 'sim time = ')
             logger.debug(self._sim_time)
             logger.debug(' haines time:')
             logger.debug( self._cluster_data['haines_time'])
             # The code really messes up when the [0] isn't in the next line. sim_time seems to just point to a place on the haines time array
-            self._sim_time = min(self._cluster_data['haines_time'])
+            self._sim_time = np.min(self._cluster_data['haines_time'])
             logger.debug( 'sim time after update= ')
             logger.debug(self._sim_time)
         else:
@@ -371,7 +371,7 @@ class InvasionPercolationTimed(GenericAlgorithm):
 
         # Mark throat as invaded
         self['throat.inv_seq'][tinvade] = self._tseq
-        self['throat.inv_pres'][tinvade] = max(max(self['throat.inv_pres']),self['throat.inv_Pc'][tinvade])
+        self['throat.inv_pres'][tinvade] = np.max((np.max(self['throat.inv_pres']),self['throat.inv_Pc'][tinvade]))
         if self._timing:
             self['throat.inv_time'][tinvade] = self._sim_time
             # update self._cluster_data.['pore_volume']
@@ -379,7 +379,7 @@ class InvasionPercolationTimed(GenericAlgorithm):
             # Remove throat's contribution to the vol_coef
             self._cluster_data['vol_coef'][self._current_cluster-1] = self._cluster_data['vol_coef'][self._current_cluster-1] - self._Tvol_coef[tinvade]
         # Mark pore as invaded
-        Pores = self._net.find_connected_pores(tinvade)
+        Pores = self._net.find_connected_pores(tinvade,flatten=True)
         # If both pores are already invaded:
         if np.in1d(Pores,np.nonzero(self['pore.cluster_final'])[0]).all():
             self._NewPore = -1
@@ -388,12 +388,12 @@ class InvasionPercolationTimed(GenericAlgorithm):
             clusters = self._cluster_data['transform'][self['pore.cluster_final'][Pores]-1]
             logger.debug('clusters = ')
             logger.debug(clusters)
-            self._current_cluster = min(clusters)
+            self._current_cluster = np.min(clusters)
             self['throat.cluster_final'][tinvade] = self._current_cluster
             # if pores are from 2 different clusters:
             if self['pore.cluster_final'][Pores[0]]!=self['pore.cluster_final'][Pores[1]] :
                 # find name of larger cluster number
-                maxCluster = max(clusters)
+                maxCluster = np.max(clusters)
                 curCluster = self._current_cluster
                 if emptyCluster == maxCluster:
                     fullCluster = curCluster
@@ -456,7 +456,7 @@ class InvasionPercolationTimed(GenericAlgorithm):
             if self._timing:
                 self['pore.inv_time'][self._NewPore] = self._sim_time
             self['pore.inv_seq'][self._NewPore] = self._tseq
-            self['pore.inv_pres'][self._NewPore] = max(self['throat.inv_pres'])
+            self['pore.inv_pres'][self._NewPore] = np.max(self['throat.inv_pres'])
             if self._timing:
                 # update self._cluster_data.['pore_volume']
                 self._cluster_data['pore_volume'][self._current_cluster-1] += self._pore_volumes[self._NewPore]
@@ -642,18 +642,18 @@ class InvasionPercolationTimed(GenericAlgorithm):
             if IPsat is not None:
                 sat_pores = self['pore.inv_sat']<=IPsat
                 sat_throats = self['throat.inv_sat']<=IPsat
-                if sum(sat_pores) == 0:
+                if np.sum(sat_pores) == 0:
                     IPseq = 0
                 else:
-                    IPseq = max([max(self['throat.inv_seq'][sat_throats]),max(self['pore.inv_seq'][sat_pores])])
+                    IPseq = np.max([np.max(self['throat.inv_seq'][sat_throats]),np.max(self['pore.inv_seq'][sat_pores])])
             else:
                 if IPpres != None:
                     sat_pores = self['pore.inv_pres']<=IPpres
                     sat_throats = self['throat.inv_pres']<=IPpres
-                    if sum(sat_pores) == 0:
+                    if np.sum(sat_pores) == 0:
                         IPseq = 0
                     else:
-                        IPseq = max([max(self['throat.inv_seq'][sat_throats]),max(self['pore.inv_seq'][sat_pores])])
+                        IPseq = np.max([np.max(self['throat.inv_seq'][sat_throats]),np.max(self['pore.inv_seq'][sat_pores])])
                 else:
                     IPseq = self._tseq
 
@@ -666,7 +666,7 @@ class InvasionPercolationTimed(GenericAlgorithm):
             inv_throats = (self['throat.inv_seq']>0)&(self['throat.inv_seq']<=IPseq)
             self._phase['throat.'+occupancy][inv_throats] = 1.
             self['throat.invaded'] = inv_throats
-            self.sat = max(self['throat.inv_sat'][inv_throats])
+            self.sat = np.max(self['throat.inv_sat'][inv_throats])
 
         except:
             print('Something bad happened while trying to update phase',self._phase.name)
diff --git a/OpenPNM/Algorithms/__OrdinaryPercolation__.py b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
index 846575f50..88bd18afe 100644
--- a/OpenPNM/Algorithms/__OrdinaryPercolation__.py
+++ b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
@@ -213,7 +213,7 @@ class OrdinaryPercolation(GenericAlgorithm):
         tmask = (pmask[temp[:, 0]] + pmask[temp[:, 1]]) * Tinvaded
         self._t_inv[(self._t_inv == sp.inf) * (tmask)] = inv_val
 
-    def evaluate_trapping(self, outlets):
+    def evaluate_trapping_tomt(self, outlets):
         r"""
         Finds trapped pores and throats after a full ordinary
         percolation drainage has been run
@@ -261,6 +261,67 @@ class OrdinaryPercolation(GenericAlgorithm):
         self['pore.inv_Pc'] = self._p_inv
         self['throat.inv_Pc'] = self._t_inv
 
+    def evaluate_trapping(self, outlets):
+        r"""
+        Finds trapped pores and throats after a full ordinary
+        percolation drainage has been run
+        Parameters
+        ----------
+        outlets : array_like
+            A list of pores that define the wetting phase outlets.
+            Disconnection from these outlets results in trapping.
+        """
+        self._p_trap = sp.zeros_like(self._p_inv, dtype=float)
+        self._t_trap = sp.zeros_like(self._t_inv, dtype=float)
+        try:
+            inv_points = sp.unique(self._p_inv)  # Get points used in OP
+        except:
+            logger.error('Orindary percolation has not been run!')
+            raise Exception('Aborting algorithm')
+        tind = self._net.throats()
+        conns = self._net.find_connected_pores(tind)
+        for inv_val in inv_points[0:-1]:
+            # Find clusters of defender pores
+            Pinvaded = self._p_inv <= inv_val
+            Cstate = sp.sum(Pinvaded[conns], axis=1)
+            Tinvaded = self._t_inv <= inv_val
+            Cstate = Cstate + Tinvaded
+            # 0 = all open, 1=1 pore filled
+            # 2=2 pores filled
+            # 3=2 pores + 1 throat filled
+            clusters = self._net.find_clusters(Cstate == 0)
+            # Clean up clusters (invaded = -1, defended >=0)
+            clusters = clusters*(~Pinvaded) - (Pinvaded)
+            # Identify clusters connected to outlet sites
+            out_clusters = sp.unique(clusters[outlets])
+            trapped_pores = ~sp.in1d(clusters, out_clusters)
+            self._p_trap[(self._p_trap == 0)[trapped_pores]] = inv_val
+            trapped_throats = self._net.find_neighbor_throats(trapped_pores,
+                                                              mode='intersection')
+            if len(trapped_throats) > 0:
+                self._t_trap[(self._t_trap == 0)[trapped_throats]] = inv_val
+        self._p_inv[self._p_trap > 0] = sp.inf
+        self._t_inv[self._t_trap > 0] = sp.inf
+        self['pore.inv_Pc'] = self._p_inv
+        self['throat.inv_Pc'] = self._t_inv
+
+    def evaluate_late_pore_filling(self, Pc, Swp_init=0.75, eta=3.0,
+                                   wetting_phase=False):
+        r"""
+        Compute the volume fraction of the phase in each pore given an initial
+        wetting phase fraction (Swp_init) and a growth exponent (eta)
+        returns the fraction of the pore volume occupied by wetting or
+        non-wetting phase.
+        Assumes Non-wetting phase displaces wetting phase
+        """
+        Swp = Swp_init*(self._p_inv/Pc)**eta
+        Swp[self._p_inv > Pc] = 1.0
+        Snwp = 1-Swp
+        if wetting_phase:
+            return Swp
+        else:
+            return Snwp
+
     def return_results(self, Pc=0, seq=None, sat=None, occupancy='occupancy'):
         r"""
         Updates the occupancy status of invading and defending phases
diff --git a/OpenPNM/Algorithms/__Tortuosity__.py b/OpenPNM/Algorithms/__Tortuosity__.py
index 1690ac4e5..754393699 100644
--- a/OpenPNM/Algorithms/__Tortuosity__.py
+++ b/OpenPNM/Algorithms/__Tortuosity__.py
@@ -16,6 +16,7 @@ import scipy.sparse.csgraph as spgr
 from OpenPNM.Algorithms import GenericAlgorithm
 import OpenPNM.Network
 from OpenPNM.Base import logging
+import OpenPNM.Utilities.misc as misc
 logger = logging.getLogger(__name__)
 
 
@@ -45,16 +46,17 @@ class Tortuosity(GenericAlgorithm):
               will require: ', t_est, ' seconds')
         return
 
-    def run(self, phase=None):
+    def run(self, phase=None, throats=None):
         logger.warning('This algorithm can take some time...')
-        graph = self._net.create_adjacency_matrix(data=self._net['throat.length'],
+        conduit_lengths = sp.sum(misc.conduit_lengths(network=self._net,
+                                 mode='centroid'), axis=1)
+        graph = self._net.create_adjacency_matrix(data=conduit_lengths,
                                                   sprsfmt='csr')
 
         if phase is not None:
             self._phase = phase
             if 'throat.occupancy' in self._phase.props():
-                temp = self._net['throat.length'] * \
-                    (self._phase['throat.occupancy'] == 1)
+                temp = conduit_lengths*(self._phase['throat.occupancy'] == 1)
                 graph = self._net.create_adjacency_matrix(data=temp,
                                                           sprsfmt='csr',
                                                           prop='temp')
diff --git a/OpenPNM/Geometry/__GenericGeometry__.py b/OpenPNM/Geometry/__GenericGeometry__.py
index 96c94be29..cf40c2a15 100644
--- a/OpenPNM/Geometry/__GenericGeometry__.py
+++ b/OpenPNM/Geometry/__GenericGeometry__.py
@@ -92,9 +92,9 @@ class GenericGeometry(Core):
             from the object.  Options are 'add' (default) or 'remove'.
 
         """
-        if len(pores) > 0:
+        if sp.size(pores) > 0:
             pores = sp.array(pores, ndmin=1)
             self._set_locations(element='pore', locations=pores, mode=mode)
-        if len(throats) > 0:
+        if sp.size(throats) > 0:
             throats = sp.array(throats, ndmin=1)
             self._set_locations(element='throat', locations=throats, mode=mode)
diff --git a/OpenPNM/Geometry/__Toray120__.py b/OpenPNM/Geometry/__Toray120__.py
new file mode 100644
index 000000000..b18758b53
--- /dev/null
+++ b/OpenPNM/Geometry/__Toray120__.py
@@ -0,0 +1,64 @@
+"""
+module __Toray120__: Subclass of GenericGeometry for a standard Toray TGPH120
+gas diffusion layer.s
+===============================================================================
+
+.. warning:: The classes of this module should be loaded through the
+'Geometry.__init__.py' file.
+
+"""
+
+from OpenPNM.Geometry import models as gm
+from OpenPNM.Geometry import GenericGeometry
+
+
+class Toray120(GenericGeometry):
+    r"""
+    Toray120 subclass of GenericGeometry
+    """
+
+    def __init__(self, **kwargs):
+        r"""
+        Initialize
+        """
+        super(Toray120, self).__init__(**kwargs)
+        self._generate()
+
+    def _generate(self):
+        r'''
+        '''
+        self.models.add(propname='pore.seed',
+                        model=gm.pore_misc.random,
+                        num_range=[0, 0.95],
+                        seed=self._seed,
+                        regen_mode='constant')
+        self.models.add(propname='throat.seed',
+                        model=gm.throat_misc.neighbor,
+                        pore_prop='pore.seed',
+                        mode='min')
+        self.models.add(propname='pore.diameter',
+                        model=gm.pore_diameter.sphere,
+                        psd_name='weibull_min',
+                        psd_shape=1.77,
+                        psd_loc=1.0e-5,
+                        psd_scale=5.0e-5,
+                        psd_offset=1.0e-6)
+        self.models.add(propname='pore.area',
+                        model=gm.pore_area.spherical)
+        self.models.add(propname='pore.volume',
+                        model=gm.pore_volume.sphere)
+        self.models.add(propname='throat.diameter',
+                        model=gm.throat_diameter.cylinder,
+                        tsd_name='weibull_min',
+                        tsd_shape=1.77,
+                        tsd_loc=5.0e-6,
+                        tsd_scale=2.5e-5,
+                        tsd_offset=1.0e-6)
+        self.models.add(propname='throat.length',
+                        model=gm.throat_length.straight)
+        self.models.add(propname='throat.volume',
+                        model=gm.throat_volume.cylinder)
+        self.models.add(propname='throat.area',
+                        model=gm.throat_area.cylinder)
+        self.models.add(propname='throat.surface_area',
+                        model=gm.throat_surface_area.cylinder)
diff --git a/OpenPNM/Geometry/__Voronoi__.py b/OpenPNM/Geometry/__Voronoi__.py
index c428ace1c..2d422cf25 100644
--- a/OpenPNM/Geometry/__Voronoi__.py
+++ b/OpenPNM/Geometry/__Voronoi__.py
@@ -9,8 +9,15 @@ Voronoi Diagram Used with Delaunay Network but could work for others (not tested
 """
 
 import scipy as sp
+import numpy as np
 from OpenPNM.Geometry import models as gm
+import OpenPNM.Utilities.vertexops as vo
 from OpenPNM.Geometry import GenericGeometry
+from OpenPNM.Base import logging
+import matplotlib.pyplot as plt
+from scipy.io import savemat
+from scipy.stats import itemfreq
+logger = logging.getLogger(__name__)
 
 
 class Voronoi(GenericGeometry):
@@ -19,15 +26,33 @@ class Voronoi(GenericGeometry):
 
     Parameters
     ----------
+    name : string
+        A unique name for the network
 
+    fibre_rad: float
+        Fibre radius to apply to Voronoi edges when calculating pore and throat
+        sizes
+
+    voxel_vol : boolean
+        Determines whether to calculate pore volumes by creating a voxel image
+        or to use the offset vertices of the throats. Voxel method is slower
+        and may run into memory issues but is more accurate and allows manipulation
+        of the image. N.B. many of the class methods are dependent on the voxel
+        image.
     """
 
-    def __init__(self, fibre_rad=3e-06, load_gen='gen', **kwargs):
+    def __init__(self, fibre_rad=3e-06, voxel_vol=True, **kwargs):
         super().__init__(**kwargs)
-        if load_gen == 'gen':
-            self._generate(fibre_rad)
+        self._fibre_rad = fibre_rad
+        self._voxel_vol = voxel_vol
+        try:
+            self._vox_len = kwargs['vox_len']
+        except:
+            self._vox_len = 1e-6
+        self._generate()
 
-    def _generate(self, fibre_rad):
+    def _generate(self):
+        # Set all the required models
         self.models.add(propname='pore.vertices',
                         model=gm.pore_vertices.voronoi)
         self.models.add(propname='throat.vertices',
@@ -36,11 +61,30 @@ class Voronoi(GenericGeometry):
                         model=gm.throat_normal.voronoi)
         self.models.add(propname='throat.offset_vertices',
                         model=gm.throat_offset_vertices.distance_transform,
-                        offset=fibre_rad,
+                        offset=self._fibre_rad,
                         set_dependent=True)
-        # Remove throats that are fully occluded and pores with no connections
-        self._net.trim_occluded_throats()
 
+        " Get rid of occluded throats in this Geometry"
+        zero_ts = self.throats()[self['throat.area'] == 0.0]
+        net_zero_ts = self.map_throats(target=self._net,
+                                       throats=zero_ts,
+                                       return_mapping=False)
+        self._net.trim(throats=net_zero_ts)
+        " Check network health and trim isolated pores "
+        h = self._net.check_network_health()
+        if np.size(h['trim_pores']) > 0:
+            self._net.trim(pores=h['trim_pores'])
+
+        if self._voxel_vol:
+            self.models.add(propname='pore.volume',
+                            model=gm.pore_volume.in_hull_volume,
+                            fibre_rad=self._fibre_rad,
+                            vox_len=self._vox_len)
+        else:
+            self.models.add(propname='pore.volume',
+                            model=gm.pore_volume.voronoi)
+        self.models.add(propname='throat.shape_factor',
+                        model=gm.throat_shape_factor.compactness)
         self.models.add(propname='pore.seed',
                         model=gm.pore_misc.random,
                         seed=self._seed)
@@ -48,26 +92,277 @@ class Voronoi(GenericGeometry):
                         model=gm.throat_misc.neighbor,
                         pore_prop='pore.seed',
                         mode='min')
-        self.models.add(propname='pore.volume',
-                        model=gm.pore_volume.voronoi)
+        self.models.add(propname='pore.centroid',
+                        model=gm.pore_centroid.voronoi)
         self.models.add(propname='pore.diameter',
-                        model=gm.pore_diameter.voronoi)
+                        model=gm.pore_diameter.equivalent_sphere)
+        self.models.add(propname='pore.indiameter',
+                        model=gm.pore_diameter.centroids)
         self.models.add(propname='pore.area',
                         model=gm.pore_area.spherical)
         self.models.add(propname='throat.diameter',
-                        model=gm.throat_diameter.voronoi)
+                        model=gm.throat_diameter.equivalent_circle)
         self.models.add(propname='throat.length',
                         model=gm.throat_length.constant,
-                        const=fibre_rad*2)
+                        const=self._fibre_rad*2)
         self.models.add(propname='throat.volume',
                         model=gm.throat_volume.extrusion)
         self.models.add(propname='throat.surface_area',
                         model=gm.throat_surface_area.extrusion)
+        self.models.add(propname='throat.c2c',
+                        model=gm.throat_length.c2c)
+
+    def make_fibre_image(self, fibre_rad=None, vox_len=1e-6):
+        r"""
+        If the voronoi voxel method was implemented to calculate pore volumes
+        an image of the fibre space has already been calculated and stored on
+        the geometry. If not generate it
+
+        Parameters
+        ----------
+        fibre_rad : float
+        Fibre radius to apply to Voronoi edges when calculating pore and throat
+        sizes
+
+        vox_len : float
+        Length of voxel edge when dividing domain
+        """
+
+        if hasattr(self, '_fibre_image'):
+            logger.info('fibre image already created')
+            return
+        else:
+            if fibre_rad is None:
+                fibre_rad = self._fibre_rad
+            fibre_rad /= vox_len
+            self._fibre_image = gm.pore_volume._get_fibre_image(self._net,
+                                                                self.pores(),
+                                                                vox_len,
+                                                                fibre_rad)
+
+    def export_fibre_image(self, mat_file='OpenPNMFibres'):
+        r"""
+        If the voronoi voxel method was implemented to calculate pore volumes
+        an image of the fibre space has already been calculated and stored on
+        the geometry.
+
+        Parameters
+        ----------
+        mat_file : string
+        Filename of Matlab file to save fibre image
+        """
+        if hasattr(self, '_fibre_image') == False:
+            logger.warning('This method only works when a fibre image exists, ' +
+                           'please run make_fibre_image')
+            return
+        matlab_dict = {"fibres": self._fibre_image}
+        savemat(mat_file, matlab_dict, format='5', long_field_names=True)
+
+    def get_fibre_slice(self, plane=None, index=None):
+        r"""
+        Plot an image of a slice through the fibre image
+        plane contains percentage values of the length of the image in each axis
+
+        Parameters
+        ----------
+        plane : array_like
+        List of 3 values, [x,y,z], 2 must be zero and the other must be between
+        zero and one representing the fraction of the domain to slice along
+        the non-zero axis
+
+        index : array_like
+        similar to plane but instead of the fraction an index of the image is used
+        """
+        if hasattr(self, '_fibre_image') == False:
+            logger.warning('This method only works when a fibre image exists, ' +
+                           'please run make_fibre_image')
+            return None
+        if plane is None and index is None:
+            logger.warning('Please provide either a plane array or index array')
+            return None
+        if self._fibre_image is None:
+            self.make_fibre_image()
+
+        if plane is not None:
+            if 'array' not in plane.__class__.__name__:
+                plane = sp.asarray(plane)
+            if sp.sum(plane == 0) != 2:
+                logger.warning('Plane argument must have two zero valued ' +
+                               'elements to produce a planar slice')
+                return None
+            l = sp.asarray(sp.shape(self._fibre_image))
+            s = sp.around(plane*l).astype(int)
+        elif index is not None:
+            if 'array' not in index.__class__.__name__:
+                index = sp.asarray(index)
+            if sp.sum(index == 0) != 2:
+                logger.warning('Index argument must have two zero valued ' +
+                               'elements to produce a planar slice')
+                return None
+            if 'int' not in str(index.dtype):
+                index = sp.around(index).astype(int)
+            s = index
+
+        if s[0] != 0:
+            slice_image = self._fibre_image[s[0], :, :]
+        elif s[1] != 0:
+            slice_image = self._fibre_image[:, s[1], :]
+        else:
+            slice_image = self._fibre_image[:, :, s[2]]
 
-if __name__ == '__main__':
-    import OpenPNM
-    pn = OpenPNM.Network.Delaunay(name='test_net')
-    pn.generate(num_pores=100, domain_size=[0.0001, 0.0001, 0.0001])
-    pn.add_boundaries()
-    test = OpenPNM.Geometry.Voronoi(pores=pn.Ps, throats=pn.Ts, network=pn)
-    pn.regenerate_geometries()
+        return slice_image
+
+    def plot_fibre_slice(self, plane=None, index=None, fig=None):
+        r"""
+        Plot one slice from the fibre image
+
+        Parameters
+        ----------
+        plane : array_like
+        List of 3 values, [x,y,z], 2 must be zero and the other must be between
+        zero and one representing the fraction of the domain to slice along
+        the non-zero axis
+
+        index : array_like
+        similar to plane but instead of the fraction an index of the image is used
+        """
+        if hasattr(self, '_fibre_image') == False:
+            logger.warning('This method only works when a fibre image exists, ' +
+                           'please run make_fibre_image')
+            return
+        slice_image = self.get_fibre_slice(plane, index)
+        if slice_image is not None:
+            if fig is None:
+                plt.figure()
+            plt.imshow(slice_image.T, cmap='Greys', origin='lower',
+                       interpolation='nearest')
+
+        return fig
+
+    def plot_porosity_profile(self, fig=None):
+        r"""
+        Return a porosity profile in all orthogonal directions by summing
+        the voxel volumes in consectutive slices.
+        """
+        if hasattr(self, '_fibre_image') == False:
+            logger.warning('This method only works when a fibre image exists, ' +
+                           'please run make_fibre_image')
+            return
+        if self._fibre_image is None:
+            self.make_fibre_image()
+        l = sp.asarray(sp.shape(self._fibre_image))
+        px = sp.zeros(l[0])
+        py = sp.zeros(l[1])
+        pz = sp.zeros(l[2])
+
+        for x in sp.arange(l[0]):
+            px[x] = sp.sum(self._fibre_image[x, :, :])
+            px[x] /= sp.size(self._fibre_image[x, :, :])
+        for y in sp.arange(l[1]):
+            py[y] = sp.sum(self._fibre_image[:, y, :])
+            py[y] /= sp.size(self._fibre_image[:, y, :])
+        for z in sp.arange(l[2]):
+            pz[z] = sp.sum(self._fibre_image[:, :, z])
+            pz[z] /= sp.size(self._fibre_image[:, :, z])
+
+        if fig is None:
+            fig = plt.figure()
+        ax = fig.gca()
+        plots = []
+        plots.append(plt.plot(sp.arange(l[0])/l[0], px, 'r', label='x'))
+        plots.append(plt.plot(sp.arange(l[1])/l[1], py, 'g', label='y'))
+        plots.append(plt.plot(sp.arange(l[2])/l[2], pz, 'b', label='z'))
+        plt.xlabel('Normalized Distance')
+        plt.ylabel('Porosity')
+        handles, labels = ax.get_legend_handles_labels()
+        ax.legend(handles, labels, loc=1)
+        plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
+        return fig
+
+    def compress_geometry(self, factor=None, preserve_fibres=False):
+        r"""
+        Adjust the vertices and recalculate geometry. Save fibre voxels before
+        and after then put them back into the image to preserve fibre volume.
+        Shape will not be conserved. Also make adjustments to the pore and throat
+        properties given an approximate volume change from adding the missing fibre
+        voxels back in
+
+        Parameters
+        ----------
+        factor : array_like
+        List of 3 values, [x,y,z], 2 must be one and the other must be between
+        zero and one representing the fraction of the domain height to compress
+        to.
+
+        preserve_fibres : boolean
+        If the fibre image has been generated and used to calculate pore volumes
+        then preserve fibre volume artificially by adjusting pore and throat sizes
+        """
+        if factor is None:
+            logger.warning('Please supply a compression factor in the form ' +
+                           '[1,1,CR], with CR < 1')
+            return
+        if sp.any(sp.asarray(factor) > 1):
+            logger.warning('The supplied compression factor is greater than 1, ' +
+                           'the method is not tested for stretching')
+            return
+        # uncompressed number of fibre voxels in each pore
+        fvu = self["pore.fibre_voxels"]
+        r1 = self["pore.diameter"]/2
+        # Most likely boundary pores - prevents divide by zero (vol change zero)
+        r1[r1 == 0] = 1
+        vo.scale(network=self._net, scale_factor=factor, preserve_vol=False)
+        self.models.regenerate()
+
+        if preserve_fibres and self._voxel_vol:
+            # compressed number of fibre voxels in each pore
+            fvc = self["pore.fibre_voxels"]
+            # amount to adjust pore volumes by
+            # (based on 1 micron cubed voxels for volume calc)
+            vol_diff = (fvu-fvc)*1e-18
+            # don't account for positive volume changes
+            vol_diff[vol_diff < 0] = 0
+            pv1 = self["pore.volume"].copy()
+            self["pore.volume"] -= vol_diff
+            self["pore.volume"][self["pore.volume"] < 0.0] = 0.0
+            pv2 = self["pore.volume"].copy()
+            "Now need to adjust the pore diameters"
+            from scipy.special import cbrt
+            rdiff = cbrt(3*np.abs(vol_diff)/(4*sp.pi))
+
+            self["pore.diameter"] -= 2*rdiff*sp.sign(vol_diff)
+            "Now as a crude approximation adjust all the throat areas and diameters"
+            "by the same ratio as the increase in a spherical pore surface area"
+            spd = np.ones(len(fvu))
+            spd[fvu > 0] = (pv2[fvu > 0]/pv1[fvu > 0])**(2/3)
+            spd[spd > 1.0] = 1.0
+            tconns = self._net["throat.conns"][self.map_throats(self._net,
+                                                                self.throats())]
+            # Need to work out the average volume change for the two pores
+            # connected by each throat Boundary pores will be connected to a
+            # throat outside this geometry if there are multiple geoms so get mapping
+            mapping = self._net.map_pores(self, self._net.pores(),
+                                          return_mapping=True)
+            source = list(mapping['source'])
+            target = list(mapping['target'])
+            ta_diff_avg = np.ones(len(tconns))
+            for i in np.arange(len(tconns)):
+                np1, np2 = tconns[i]
+                if np1 in source and np2 in source:
+                    gp1 = target[source.index(np1)]
+                    gp2 = target[source.index(np2)]
+                    ta_diff_avg[i] = (spd[gp1] + spd[gp2]) / 2
+                elif np1 in source:
+                    gp1 = target[source.index(np1)]
+                    ta_diff_avg[i] = spd[gp1]
+                elif np2 in source:
+                    gp2 = target[source.index(np2)]
+                    ta_diff_avg[i] = spd[gp2]
+            self["throat.area"] *= ta_diff_avg
+            self["throat.area"][self["throat.area"] < 0] = 0
+            self["throat.diameter"] = 2*sp.sqrt(self["throat.area"]/sp.pi)
+            self["throat.indiameter"] *= sp.sqrt(ta_diff_avg)
+        else:
+            logger.warning('Fibre volume is not be conserved under compression')
+        # Remove pores with zero throats
+        self._net.trim_occluded_throats()
diff --git a/OpenPNM/Geometry/__init__.py b/OpenPNM/Geometry/__init__.py
index f04523c86..7b41f7ba5 100644
--- a/OpenPNM/Geometry/__init__.py
+++ b/OpenPNM/Geometry/__init__.py
@@ -41,5 +41,6 @@ from .__TestGeometry__ import TestGeometry
 from .__Boundary__ import Boundary
 from .__Voronoi__ import Voronoi
 from .__Toray090__ import Toray090
+from .__Toray120__ import Toray120
 from .__SGL10__ import SGL10
 from . import models
diff --git a/OpenPNM/Geometry/models/pore_centroid.py b/OpenPNM/Geometry/models/pore_centroid.py
index 0fab83a14..fabed9e24 100644
--- a/OpenPNM/Geometry/models/pore_centroid.py
+++ b/OpenPNM/Geometry/models/pore_centroid.py
@@ -1,37 +1,29 @@
 r"""
 ===============================================================================
-pore_centroid -- N.B Neither method is accurate but is quicker than that
-                 implemented in pore.volume
+pore_centroid
 ===============================================================================
 
 """
 import scipy as _sp
 
 
-def voronoi(geometry, pore_vertices='pore.vertices', **kwargs):
-    r"""
-    Calculate the centroid of the pore from the voronoi vertices - C.O.M
-    """
-    verts = geometry[pore_vertices]
-    value = _sp.ndarray([len(verts), 3])
-    for i, vert in enumerate(verts):
-        value[i] = _sp.array([vert[:, 0].mean(),
-                              vert[:, 1].mean(),
-                              vert[:, 2].mean()])
-    return value
-
-
-def voronoi2(geometry, vertices='throat.centroid', **kwargs):
+def voronoi(network, geometry, vertices='throat.centroid', **kwargs):
     r"""
     Calculate the centroid from the mean of the throat centroids
     """
     value = _sp.ndarray([geometry.num_pores(), 3])
-    pore_map = geometry.map_pores(geometry.pores(), geometry._net)
-    for geom_pore, net_pore in pore_map:
-        net_throats = geometry._net.find_neighbor_throats(net_pore)
-        geom_throats = geometry._net.map_throats(net_throats, geometry)[:, 1]
+    pore_map = geometry.map_pores(target=network,
+                                  pores=geometry.pores(),
+                                  return_mapping=True)
+    for i, net_pore in enumerate(pore_map['target']):
+        geom_pore = pore_map['source'][i]
+        net_throats = network.find_neighbor_throats(net_pore)
+        geom_throats = network.map_throats(target=geometry,
+                                           throats=net_throats,
+                                           return_mapping=False)
         verts = geometry[vertices][geom_throats]
-        value[geom_pore] = _sp.array([verts[:, 0].mean(),
-                                      verts[:, 1].mean(),
-                                      verts[:, 2].mean()])
+        " Ignore all zero centroids "
+        verts = verts[~_sp.all(verts == 0, axis=1)]
+        value[geom_pore] = _sp.mean(verts, axis=0)
+
     return value
diff --git a/OpenPNM/Geometry/models/pore_diameter.py b/OpenPNM/Geometry/models/pore_diameter.py
index de8a77619..19bc9b00a 100644
--- a/OpenPNM/Geometry/models/pore_diameter.py
+++ b/OpenPNM/Geometry/models/pore_diameter.py
@@ -86,17 +86,6 @@ def equivalent_cube(geometry, pore_volume='pore.volume', **kwargs):
     return value
 
 
-def voronoi(geometry, pore_volume='pore.volume', **kwargs):
-    r"""
-    Calculate pore diameter from equivalent sphere - volumes must be calculated
-    first.
-    """
-    from scipy.special import cbrt
-    pore_vols = geometry[pore_volume]
-    value = cbrt(6*pore_vols/_sp.pi)
-    return value
-
-
 def centroids(network, geometry, **kwargs):
     r"""
     Calculate the diameter representing an inclosed sphere. The maximum is very
@@ -106,87 +95,19 @@ def centroids(network, geometry, **kwargs):
     """
     Np = geometry.num_pores()
     value = _sp.zeros(Np)
-    pore_map = geometry.map_pores(geometry.pores(), geometry._net)
-    for geom_pore, net_pore in pore_map:
+    pore_map = geometry.map_pores(target=network,
+                                  pores=geometry.pores(),
+                                  return_mapping=True)
+
+    for i, net_pore in enumerate(pore_map['target']):
+        geom_pore = pore_map['source'][i]
         net_throats = geometry._net.find_neighbor_throats(net_pore)
-        geom_throats = geometry._net.map_throats(net_throats, geometry)[:, 1]
+        geom_throats = geometry._net.map_throats(target=geometry,
+                                                 throats=net_throats,
+                                                 return_mapping=False)
         tcs = geometry["throat.centroid"][geom_throats]
         pc = geometry["pore.centroid"][geom_pore]
-        value[geom_pore] = _sp.mean(_sp.sqrt(((tcs - pc) * (tcs - pc))[:, 0] +
+        value[geom_pore] = _sp.mean(_sp.sqrt(((tcs-pc)*(tcs-pc))[:, 0] +
                                              ((tcs-pc)*(tcs-pc))[:, 1] +
                                              ((tcs-pc)*(tcs-pc))[:, 2]))*2
     return value
-
-
-def insphere(network, geometry, **kwargs):
-    r"""
-    Calculate the diameter of the insphere using linear programming
-    """
-    import warnings
-    try:
-        import pulp as pu
-        Np = geometry.num_pores()
-        value = _sp.zeros(Np)
-        pore_map = geometry.map_pores(geometry.pores(), geometry._net)
-        for geom_pore, net_pore in pore_map:
-            net_throats = geometry._net.find_neighbor_throats(net_pore)
-            geom_throats = geometry._net.map_throats(net_throats, geometry)[:, 1]
-            verts = geometry['throat.offset_vertices'][geom_throats]
-            if len(verts) > 1:
-                try:
-                    pts = np.vstack((i for i in verts if len(i) > 0))
-                except ValueError:
-                    pts = []
-                if len(pts) > 4:
-                    "Work out central point to use as initial guess"
-                    c0 = np.mean(pts, axis=0)
-                    "Compute convex hull to find points lying on the hull in order"
-                    hull = ConvexHull(pts, qhull_options='QJ Pp')
-                    "For each simplex making up the hull collect the end points"
-                    A = pts[hull.simplices[:, 0]]
-                    B = pts[hull.simplices[:, 1]]
-                    C = pts[hull.simplices[:, 2]]
-                    "Normal of the simplices"
-                    N = np.cross((B-A), (C-A), axis=1)
-                    "Normalize the normal vector"
-                    L = np.linalg.norm(N, axis=1)
-                    F = np.vstack((L, L, L)).T
-                    N /= F
-                    "If normals point out of hull change sign to point in"
-                    pointing_out = (np.sum((A-c0)*N, axis=1) > 0)
-                    N[pointing_out] *= -1
-                    "Define Linear Program Variables"
-                    "The centre of the incircle adjustment"
-                    cx = pu.LpVariable('cx', None, None, pu.LpContinuous)
-                    cy = pu.LpVariable('cy', None, None, pu.LpContinuous)
-                    cz = pu.LpVariable('cy', None, None, pu.LpContinuous)
-                    "Radius of the incircle"
-                    R = pu.LpVariable('R', 0, None, pu.LpContinuous)
-                    # Slack variables for shortest distance between centre
-                    # and simplices
-                    S = pu.LpVariable.dict('SlackVariable', range(len(A)), 0,
-                                           None, pu.LpContinuous)
-                    "Set up LP problem"
-                    prob = pu.LpProblem('FindInRadius', pu.LpMaximize)
-                    "Objective Function"
-                    prob += R
-                    for i in range(len(A)):
-                        "Ni.(C-Ai)-Si = 0"
-                        prob += N[i][0]*(c0[0] + cx) + N[i][1]*(c0[1] + cy) + \
-                            N[i][2]*(c0[2] + cz) - N[i][0]*A[i][0] - \
-                            N[i][1]*A[i][1] - N[i][2]*A[i][2] - S[i] == 0
-                        "Si >= R"
-                        prob += S[i] >= R
-                    "Solve the LP"
-                    with warnings.catch_warnings():
-                        warnings.simplefilter('ignore')
-                        prob.solve()
-                    # As the radius is the objective function we can get it from the
-                    # objective or as R.value()
-                    rad = prob.objective.value()
-                    value[geom_pore] = rad*2
-
-        return value
-
-    except ImportError:
-        print('Cannot use insphere method without installing pulp package')
diff --git a/OpenPNM/Geometry/models/pore_volume.py b/OpenPNM/Geometry/models/pore_volume.py
index ae7843f8e..201f8a0e1 100644
--- a/OpenPNM/Geometry/models/pore_volume.py
+++ b/OpenPNM/Geometry/models/pore_volume.py
@@ -7,7 +7,243 @@ pore_volume --
 import scipy as _sp
 import numpy as np
 from scipy.spatial import Delaunay
+from scipy.spatial import ConvexHull
 import OpenPNM.Utilities.misc as misc
+from scipy import ndimage
+from OpenPNM.Base import logging
+logger = logging.getLogger(__name__)
+
+
+def inhull(geometry, xyz, pore, tol=1e-7):
+    r"""
+    Tests whether points lie within a convex hull or not.
+    Computes a tesselation of the hull works out the normals of the facets.
+    Then tests whether dot(x.normals) < dot(a.normals) where a is the the
+    first vertex of the facets
+    """
+    xyz = np.around(xyz, 10)
+    " Work out range to span over for pore hull "
+    xmin = xyz[:, 0].min()
+    xr = (np.ceil(xyz[:, 0].max())-np.floor(xmin)).astype(int)+1
+    ymin = xyz[:, 1].min()
+    yr = (np.ceil(xyz[:, 1].max())-np.floor(ymin)).astype(int)+1
+    zmin = xyz[:, 2].min()
+    zr = (np.ceil(xyz[:, 2].max())-np.floor(zmin)).astype(int)+1
+
+    origin = np.array([xmin, ymin, zmin])
+    " start index "
+    si = np.floor(origin).astype(int)
+    xyz -= origin
+    dom = np.zeros([xr, yr, zr], dtype=np.uint8)
+    indx, indy, indz = np.indices((xr, yr, zr))
+    " Calculate the tesselation of the points "
+    hull = ConvexHull(xyz)
+    " Assume 3d for now "
+    " Calculate normals from the vector cross product of the vectors "
+    " defined by joining points in the simplices "
+    vab = xyz[hull.simplices[:, 0]]-xyz[hull.simplices[:, 1]]
+    vac = xyz[hull.simplices[:, 0]]-xyz[hull.simplices[:, 2]]
+    nrmls = np.cross(vab, vac)
+    " Scale normal vectors to unit length "
+    nrmlen = np.sum(nrmls**2, axis=-1)**(1./2)
+    nrmls = nrmls*np.tile((1/nrmlen), (3, 1)).T
+    " Center of Mass "
+    center = np.mean(xyz, axis=0)
+    " Any point from each simplex "
+    a = xyz[hull.simplices[:, 0]]
+    " Make sure all normals point inwards "
+    dp = np.sum((np.tile(center, (len(a), 1))-a)*nrmls, axis=-1)
+    k = dp < 0
+    nrmls[k] = -nrmls[k]
+    " Now we want to test whether dot(x,N) >= dot(a,N) "
+    aN = np.sum(nrmls*a, axis=-1)
+
+    for plane_index in range(len(a)):
+        eqx = nrmls[plane_index][0]*(indx)
+        eqy = nrmls[plane_index][1]*(indy)
+        eqz = nrmls[plane_index][2]*(indz)
+        xN = eqx + eqy + eqz
+        dom[xN - aN[plane_index] >= 0-tol] += 1
+    dom[dom < len(a)] = 0
+    dom[dom == len(a)] = 1
+    ds = np.shape(dom)
+    temp_arr = np.zeros_like(geometry._hull_image, dtype=bool)
+    temp_arr[si[0]:si[0]+ds[0], si[1]:si[1]+ds[1], si[2]:si[2]+ds[2]] = dom
+    geometry._hull_image[temp_arr] = pore
+    hull_num = np.sum(dom)
+    dom = dom * geometry._fibre_image[si[0]:si[0]+ds[0], si[1]:si[1]+ds[1],
+                                      si[2]:si[2]+ds[2]]
+    pore_num = np.sum(dom)
+    fibre_num = hull_num - pore_num
+
+    del temp_arr
+    return pore_num, fibre_num
+
+
+def _voxel_centroid(image, pores=None, vox_len=1e-6):
+    r"""
+    Calculate Pore Centroid from indices of the voronoi voxel image generated
+    in _get_voxel_volume
+    """
+    if pores is None:
+        pores = np.unique(image)
+    centroids = np.zeros([len(pores), 3])
+    for i, pore in enumerate(pores):
+        px, py, pz = np.where(image == pore)
+        try:
+            " if one is empty then all will be - no pore volume for this pore "
+            if len(px) > 0:
+                cx = np.mean(px)
+                cy = np.mean(py)
+                cz = np.mean(pz)
+            centroids[i] = np.array([cx, cy, cz])*vox_len
+        except:
+            logger.warn("Some pore centroid data may be invalid, look for zeros")
+    return centroids
+
+
+def _get_vertex_range(verts):
+    "Find the extent of the vetrices"
+    [vxmin, vxmax, vymin, vymax, vzmin, vzmax] = [1e20, 0, 1e20, 0, 1e20, 0]
+    for vert in verts:
+        if np.min(vert[:, 0]) < vxmin:
+            vxmin = np.min(vert[:, 0])
+        if np.max(vert[:, 0]) > vxmax:
+            vxmax = np.max(vert[:, 0])
+        if np.min(vert[:, 1]) < vymin:
+            vymin = np.min(vert[:, 1])
+        if np.max(vert[:, 1]) > vymax:
+            vymax = np.max(vert[:, 1])
+        if np.min(vert[:, 2]) < vzmin:
+            vzmin = np.min(vert[:, 2])
+        if np.max(vert[:, 2]) > vzmax:
+            vzmax = np.max(vert[:, 2])
+    return [vxmin, vxmax, vymin, vymax, vzmin, vzmax]
+
+
+def _get_fibre_image(network, cpores, vox_len, fibre_rad):
+    r"""
+    Produce image by filling in voxels along throat edges using Bresenham line
+    Then performing distance transform on fibre voxels to erode the pore space
+    """
+
+    cthroats = network.find_neighbor_throats(pores=cpores)
+
+    "Below method copied from geometry model throat.vertices"
+    "Needed now as network may not have all throats assigned to geometry"
+    "i.e network['throat.vertices'] could return garbage"
+    verts = _sp.ndarray(network.num_throats(), dtype=object)
+    for i in range(len(verts)):
+        verts[i] = _sp.asarray(list(network["throat.vert_index"][i].values()))
+    cverts = verts[cthroats]
+    [vxmin, vxmax, vymin, vymax, vzmin, vzmax] = _get_vertex_range(cverts)
+    "Translate vertices so that minimum occurs at the origin"
+    for index in range(len(cverts)):
+        cverts[index] -= np.array([vxmin, vymin, vzmin])
+    "Find new size of image array"
+    cdomain = np.around(np.array([(vxmax-vxmin), (vymax-vymin), (vzmax-vzmin)]), 6)
+    logger.info("Creating fibre domain range: " + str(np.around(cdomain, 5)))
+    lx = np.int(np.around(cdomain[0]/vox_len)+1)
+    ly = np.int(np.around(cdomain[1]/vox_len)+1)
+    lz = np.int(np.around(cdomain[2]/vox_len)+1)
+    "Try to create all the arrays we will need at total domain size"
+    try:
+        pore_space = np.ones([lx, ly, lz], dtype=np.uint8)
+        fibre_space = np.zeros(shape=[lx, ly, lz], dtype=np.uint8)
+        dt = np.zeros([lx, ly, lz], dtype=float)
+        " Only need one chunk "
+        cx = cy = cz = 1
+        chunk_len = np.max(np.shape(pore_space))
+    except:
+        logger.info("Domain too large to fit into memory so chunking domain to" +
+                    "process image, this may take some time")
+        " Do chunking "
+        chunk_len = 100
+        if (lx > chunk_len):
+            cx = np.ceil(lx/chunk_len).astype(int)
+        else:
+            cx = 1
+        if (ly > chunk_len):
+            cy = np.ceil(ly/chunk_len).astype(int)
+        else:
+            cy = 1
+        if (lz > chunk_len):
+            cz = np.ceil(lz/chunk_len).astype(int)
+        else:
+            cz = 1
+
+    "Get image of the fibres"
+    line_points = bresenham(cverts, vox_len/2)
+    line_ints = (np.around((line_points/vox_len), 0)).astype(int)
+    for x, y, z in line_ints:
+        try:
+            pore_space[x][y][z] = 0
+        except IndexError:
+            logger.warning("Some elements in image processing are out of bounds")
+
+    num_chunks = np.int(cx*cy*cz)
+    cnum = 1
+    for ci in range(cx):
+        for cj in range(cy):
+            for ck in range(cz):
+                "Work out chunk range"
+                logger.info("Processing Fibre Chunk: "+str(cnum)+" of " +
+                            str(num_chunks))
+                cxmin = ci*chunk_len
+                cxmax = np.int(np.ceil((ci+1)*chunk_len + 5*fibre_rad))
+                cymin = cj*chunk_len
+                cymax = np.int(np.ceil((cj+1)*chunk_len + 5*fibre_rad))
+                czmin = ck*chunk_len
+                czmax = np.int(np.ceil((ck+1)*chunk_len + 5*fibre_rad))
+                "Don't overshoot"
+                if cxmax > lx:
+                    cxmax = lx
+                if cymax > ly:
+                    cymax = ly
+                if czmax > lz:
+                    czmax = lz
+                dt = ndimage.distance_transform_edt(pore_space[cxmin:cxmax,
+                                                               cymin:cymax,
+                                                               czmin:czmax])
+                fibre_space[cxmin:cxmax,
+                            cymin:cymax,
+                            czmin:czmax][dt <= fibre_rad] = 0
+                fibre_space[cxmin:cxmax,
+                            cymin:cymax,
+                            czmin:czmax][dt > fibre_rad] = 1
+                del dt
+                cnum += 1
+    del pore_space
+    return fibre_space
+
+
+def bresenham(faces, dx):
+    line_points = []
+    for face in faces:
+        " Get in hull order "
+        fx = face[:, 0]
+        fy = face[:, 1]
+        fz = face[:, 2]
+        " Find the axis with the smallest spread of points and remove it to make 2D "
+        if (np.std(fx) < np.std(fy)) and (np.std(fx) < np.std(fz)):
+            f2d = np.vstack((fy, fz)).T
+        elif (np.std(fy) < np.std(fx)) and (np.std(fy) < np.std(fz)):
+            f2d = np.vstack((fx, fz)).T
+        else:
+            f2d = np.vstack((fx, fy)).T
+        hull = ConvexHull(f2d, qhull_options='QJ Pp')
+        face = np.around(face[hull.vertices], 6)
+        for i in range(len(face)):
+            vec = face[i]-face[i-1]
+            vec_length = np.linalg.norm(vec)
+            increments = np.ceil(vec_length/dx)
+            check_p_old = np.array([-1, -1, -1])
+            for x in np.linspace(0, 1, increments):
+                check_p_new = face[i-1]+(vec*x)
+                if np.sum(check_p_new - check_p_old) != 0:
+                    line_points.append(check_p_new)
+                    check_p_old = check_p_new
+    return np.asarray(line_points)
 
 
 def sphere(geometry, pore_diameter='pore.diameter', **kwargs):
@@ -56,58 +292,60 @@ def _get_hull_volume(points):
     into triangles and working out the volume of all the pyramid elements connected
     to the volume centroid
     """
-    # Remove any duplicate points - this messes up the triangulation
+    " Remove any duplicate points - this messes up the triangulation "
     points = _sp.asarray(misc.unique_list(np.around(points, 10)))
     try:
         tri = Delaunay(points, qhull_options='QJ Pp')
     except _sp.spatial.qhull.QhullError:
         print(points)
-    # We only want points included in the convex hull to calculate the centroid
+    " We only want points included in the convex hull to calculate the centroid "
     hull_centroid = _sp.array([points[:, 0].mean(),
                                points[:, 1].mean(),
                                points[:, 2].mean()])
     hull_volume = 0.0
     pyramid_COMs = []
     for ia, ib, ic in tri.convex_hull:
-        # Points making each triangular face
-        # Collection of co-ordinates of each point in this face
+        " Points making each triangular face "
+        " Collection of co-ordinates of each point in this face "
         face_x = points[[ia, ib, ic]][:, 0]
         face_y = points[[ia, ib, ic]][:, 1]
         face_z = points[[ia, ib, ic]][:, 2]
-        # Average of each co-ordinate is the centroid of the face
+        " Average of each co-ordinate is the centroid of the face "
         face_centroid = [face_x.mean(), face_y.mean(), face_z.mean()]
         face_centroid_vector = face_centroid - hull_centroid
-        # Vectors of the sides of the face used to find normal vector and area
+        " Vectors of the sides of the face used to find normal vector and area "
         vab = points[ib] - points[ia]
         vac = points[ic] - points[ia]
         vbc = points[ic] - points[ib]
-        # As vectors are co-planar the cross-product will produce the normal
-        # vector of the face
+        " As vectors are co-planar the cross-product will produce the normal "
+        " vector of the face "
         face_normal = _sp.cross(vab, vac)
         try:
             face_unit_normal = face_normal/_sp.linalg.norm(face_normal)
         except RuntimeWarning:
             print('Pore Volume Error:' + str(vab) + ' ' + str(vac))
-        # As triangles are orientated randomly in 3D we could either transform
-        # co-ordinates to align with a plane and perform 2D operations to work out
-        # the area or we could work out the lengths of each side and use Heron's
-        # formula which is easier. Using Delaunay traingulation will always produce
-        # triangular faces but if dealing with other polygons co-ordinate transfer
-        # may be necessary
+        """
+        As triangles are orientated randomly in 3D we could either transform
+        co-ordinates to align with a plane and perform 2D operations to work out
+        the area or we could work out the lengths of each side and use Heron's
+        formula which is easier. Using Delaunay traingulation will always produce
+        triangular faces but if dealing with other polygons co-ordinate transfer
+        may be necessary
+        """
         a = _sp.linalg.norm(vab)
         b = _sp.linalg.norm(vbc)
         c = _sp.linalg.norm(vac)
-        # Semiperimeter
+        " Semiperimeter "
         s = 0.5*(a + b + c)
         face_area = _sp.sqrt(s*(s-a)*(s-b)*(s-c))
-        # Now the volume of the pyramid section defined by the 3 face points and the
-        # hull centroid can be calculated
+        " Now the volume of the pyramid section defined by the 3 face points and "
+        " the hull centroid can be calculated "
         pyramid_volume = _sp.absolute(_sp.dot(face_centroid_vector,
                                               face_unit_normal)*face_area/3)
-        # Each pyramid is summed together to calculate the total volume
+        " Each pyramid is summed together to calculate the total volume "
         hull_volume += pyramid_volume
-        # The Centre of Mass will not be the same as the geometrical centroid. A
-        # weighted adjustment can be calculated from the pyramid centroid and volume
+        " The Centre of Mass will not be the same as the geometrical centroid. "
+        " Weighted adjustment is calculated from the pyramid centroid and volume "
         vha = points[ia] - hull_centroid
         vhb = points[ib] - hull_centroid
         vhc = points[ic] - hull_centroid
@@ -153,8 +391,8 @@ def voronoi(network, geometry, **kwargs):
         elif len(geom_throats) == 1 and 'throat.centroid' in geometry.props():
                 com[i] = geometry['throat.centroid'][geom_throats]
                 volume[i] = 0
-    # Find any pores with centroids at origin and use the mean of the pore vertices
-    # instead. Not doing this messes up hydraulic conductances using centre to centre
+    " Find any pores with centroids at origin and use the mean of the pore vertices "
+    " Not doing this messes up hydraulic conductances using centre to centre "
     ps = np.where(~com.any(axis=1))[0]
     if len(ps) > 0:
         for pore in ps:
@@ -162,3 +400,41 @@ def voronoi(network, geometry, **kwargs):
     geometry['pore.centroid'] = com
 
     return volume
+
+
+def in_hull_volume(network, geometry, fibre_rad, vox_len=1e-6, **kwargs):
+    r"""
+    Work out the voxels inside the convex hull of the voronoi vertices of each pore
+    """
+    Np = network.num_pores()
+    geom_pores = geometry.map_pores(network, geometry.pores())
+    volume = _sp.zeros(Np)
+    pore_vox = _sp.zeros(Np, dtype=int)
+    fibre_vox = _sp.zeros(Np, dtype=int)
+    voxel = vox_len**3
+    try:
+        nbps = network.pores('boundary', mode='not')
+    except KeyError:
+        " Boundaries have not been generated "
+        nbps = network.pores()
+    " Voxel length "
+    fibre_rad = np.around((fibre_rad-(vox_len/2))/vox_len, 0).astype(int)
+
+    "Get the fibre image"
+    fibre_image = _get_fibre_image(network, geom_pores, vox_len, fibre_rad)
+    " Save as private variables "
+    geometry._fibre_image = fibre_image
+    hull_image = np.ones_like(fibre_image, dtype=np.uint16)*-1
+    geometry._hull_image = hull_image
+    for pore in nbps:
+        logger.info("Processing Pore: "+str(pore+1)+" of "+str(len(nbps)))
+        verts = np.asarray([i for i in network["pore.vert_index"][pore].values()])
+        verts = np.asarray(misc.unique_list(np.around(verts, 6)))
+        verts /= vox_len
+        pore_vox[pore], fibre_vox[pore] = inhull(geometry, verts, pore)
+
+    volume = pore_vox*voxel
+    geometry["pore.fibre_voxels"] = fibre_vox[geom_pores]
+    geometry["pore.pore_voxels"] = pore_vox[geom_pores]
+
+    return volume[geom_pores]
diff --git a/OpenPNM/Geometry/models/throat_area.py b/OpenPNM/Geometry/models/throat_area.py
index 21deb26d6..6dfa5b469 100644
--- a/OpenPNM/Geometry/models/throat_area.py
+++ b/OpenPNM/Geometry/models/throat_area.py
@@ -7,6 +7,7 @@ Submodule -- throat_area
 import scipy as _sp
 import OpenPNM.Utilities.transformations as tr
 import OpenPNM.Utilities.vertexops as vo
+from scipy.spatial import ConvexHull
 
 
 def cylinder(geometry, throat_diameter='throat.diameter', **kwargs):
@@ -60,6 +61,9 @@ def voronoi(geometry, **kwargs):
     for i in range(Nt):
         if len(verts[i]) > 2:
             verts_2D = tr.rotate_and_chop(verts[i], normals[i], [0, 0, 1])
+            # Get in hull order
+            hull = ConvexHull(verts_2D, qhull_options='QJ Pp')
+            verts_2D = verts_2D[hull.vertices]
             area[i] = vo.PolyArea2D(verts_2D)
         else:
             area[i] = 0.0
diff --git a/OpenPNM/Geometry/models/throat_centroid.py b/OpenPNM/Geometry/models/throat_centroid.py
index bfc4ba3ea..3cdaf0352 100644
--- a/OpenPNM/Geometry/models/throat_centroid.py
+++ b/OpenPNM/Geometry/models/throat_centroid.py
@@ -10,35 +10,13 @@ from OpenPNM.Utilities import vertexops as vo
 from scipy.spatial import ConvexHull
 
 
-def voronoi(geometry, **kwargs):
-    r"""
-    Calculate the geometrical centroid of the throat from the voronoi vertices.
-    """
-    verts = geometry['throat.vertices']
-    offset_verts = geometry['throat.offset_vertices']
-    value = _sp.ndarray([len(verts), 3])
-    for i in range(len(verts)):
-        if len(offset_verts[i]) > 2:
-            value[i] = _sp.array([offset_verts[i][:, 0].mean(),
-                                  offset_verts[i][:, 1].mean(),
-                                  offset_verts[i][:, 2].mean()])
-        elif len(verts[i]) > 2:
-            value[i] = _sp.array([verts[i][:, 0].mean(),
-                                  verts[i][:, 1].mean(),
-                                  verts[i][:, 2].mean()])
-        else:
-            value[i] = _sp.array([0, 0, 0])
-
-    return value
-
-
-def centre_of_mass(geometry, **kwargs):
+def centre_of_mass(geometry, vertices='throat.offset_vertices', **kwargs):
     r"""
     Calculate the centre of mass of the throat from the voronoi vertices.
     """
     Nt = geometry.num_throats()
     outer_verts = geometry['throat.vertices']
-    offset_verts = geometry['throat.offset_vertices']
+    offset_verts = geometry[vertices]
     normal = geometry['throat.normal']
     z_axis = [0, 0, 1]
     value = _sp.ndarray([Nt, 3])
diff --git a/OpenPNM/Geometry/models/throat_diameter.py b/OpenPNM/Geometry/models/throat_diameter.py
index 92a9f8a8b..42821a1b4 100644
--- a/OpenPNM/Geometry/models/throat_diameter.py
+++ b/OpenPNM/Geometry/models/throat_diameter.py
@@ -7,9 +7,6 @@ Submodule -- throat_diameter
 import scipy as _sp
 import numpy as np
 
-import OpenPNM.Utilities.transformations as tr
-from scipy.spatial import ConvexHull
-
 
 def cylinder(geometry, tsd_name, tsd_shape, tsd_loc, tsd_scale,
              throat_seed='throat.seed', tsd_offset=0, **kwargs):
@@ -23,81 +20,24 @@ def cylinder(geometry, tsd_name, tsd_shape, tsd_loc, tsd_scale,
     return value
 
 
-def voronoi(geometry, throat_area='throat.area', **kwargs):
+def equivalent_circle(geometry, throat_area='throat.area', **kwargs):
     r"""
-    Calculate throat diameter from analysis of Voronoi facets
-    Equivalent circular diameter from voronoi area
-    Could do better here and work out minimum diameter from verts
+    Equivalent circular diameter from throat area
     """
     areas = geometry[throat_area]
     value = 2*_sp.sqrt(areas/_sp.pi)  # 64 bit sqrt doesn't work!
     return value
 
 
-def incircle(geometry, **kwargs):
+def minpore(network, geometry, **kwargs):
     r"""
-    Calculate the incircle diameter by linear programming and the simplex search
-    algorithm using the offset vertices generated by the voronoi diagram offsetting
-    routine
+    Assign the throat diameter to be equal to the smallest connecting pore
+    diameter. If zero (in case of boundaries) take it to be the maximum of
+    the connecting pore diameters
     """
-    import warnings
-    try:
-        import pulp as pu
-        Nt = geometry.num_throats()
-        verts = geometry['throat.offset_vertices']
-        normals = geometry['throat.normal']
-        value = _sp.zeros(Nt)
-        for i in range(Nt):
-            if len(verts[i]) > 2:
-                pts = tr.rotate_and_chop(verts[i], normals[i], [0, 0, 1])
-                # Work out central point to use as initial guess
-                C = np.mean(pts, axis=0)
-                # Compute convex hull to find points lying on the hull in order
-                hull = ConvexHull(pts, qhull_options='QJ Pp')
-                # For each simplex making up the hull collect the end points
-                A = pts[hull.vertices]
-                B = pts[np.roll(hull.vertices, -1)]
-                I = np.array([[0, 1], [-1, 0]])
-                # Normal of the simplices
-                N = np.dot((B-A), I)
-                # Normalize the normal vector
-                L = np.linalg.norm(N, axis=1)
-                F = np.vstack((L, L)).T
-                N /= F
-                # Mid-points of the simplex
-                M = (B+A)/2
-                # If normals point out of hull change sign to point in
-                pointing_out = (np.sum((M-C)*N, axis=1) > 0)
-                N[pointing_out] *= -1
-                # Define Linear Program Variables
-                # The centre of the incircle adjustment
-                cx = pu.LpVariable('cx', None, None, pu.LpContinuous)
-                cy = pu.LpVariable('cy', None, None, pu.LpContinuous)
-                # Radius of the incircle
-                R = pu.LpVariable('R', 0, None, pu.LpContinuous)
-                # Slack variables for shortest distance between centre and simplices
-                S = pu.LpVariable.dict('SlackVariable', range(len(A)), 0,
-                                       None, pu.LpContinuous)
-                # Set up LP problem
-                prob = pu.LpProblem('FindInRadius', pu.LpMaximize)
-                # Objective Function
-                prob += R
-                for j in range(len(A)):
-                    # Ni.(C-Ai)-Si = 0
-                    prob += N[j][0]*(C[0]+cx) + N[j][1]*(C[1]+cy) - \
-                        N[j][0]*A[j][0] - N[j][1]*A[j][1] - S[j] == 0
-                    # Si >= R
-                    prob += S[j] >= R
-                # Solve the LP
-                with warnings.catch_warnings():
-                    warnings.simplefilter('ignore')
-                    prob.solve()
-                # As the radius is the objective function we can get it from the
-                # objective or as R.value()
-                value[i] = 2*R.value()
-            else:
-                value[i] = 0.0
-
-        return value
-    except ImportError:
-        print('Cannot use incircle method without installing pulp package')
+    gTs = geometry.throats()
+    nTs = geometry.map_throats(network, gTs)
+    pDs = network["pore.diameter"][network["throat.conns"][nTs]]
+    value = np.min(pDs, axis=1)
+    value[value == 0.0] = np.max(pDs, axis=1)[value == 0.0]
+    return value
diff --git a/OpenPNM/Geometry/models/throat_length.py b/OpenPNM/Geometry/models/throat_length.py
index 2e5bf597e..ac42214fb 100644
--- a/OpenPNM/Geometry/models/throat_length.py
+++ b/OpenPNM/Geometry/models/throat_length.py
@@ -38,7 +38,7 @@ def straight(network, geometry, pore_diameter='pore.diameter',
     return value
 
 
-def voronoi(network, geometry, **kwargs):
+def c2c(network, geometry, **kwargs):
     r"""
     Calculate the centre to centre distance from centroid of pore1 to centroid of
     throat to centroid of pore2. This is tricky as connections are defined at
@@ -61,27 +61,6 @@ def voronoi(network, geometry, **kwargs):
     return value
 
 
-def c2c(network, geometry, **kwargs):
-    r"""
-    Calculate throat length
-    """
-    # Initialize throat_property['length']
-    throats = network.throats(geometry.name)
-    pore1 = network['throat.conns'][:, 0]
-    pore2 = network['throat.conns'][:, 1]
-    C1 = network['pore.coords'][pore1]
-    C2 = network['pore.coords'][pore2]
-    E = _sp.sqrt(_sp.sum((C1-C2)**2, axis=1))  # Euclidean distance between pores
-    value = E
-    value = value[throats]
-    if _sp.any(value < 0):
-        geometry._logger.warning('Negative throat lengths are calculated. \
-                                  Arbitrary positive length assigned (1e9 meters)')
-        Ts = _sp.where(value < 0)[0]
-        value[Ts] = 1e-9
-    return value
-
-
 def constant(network, geometry, const, **kwargs):
     r"""
     Calculate throat length
diff --git a/OpenPNM/Geometry/models/throat_offset_vertices.py b/OpenPNM/Geometry/models/throat_offset_vertices.py
index 08291159d..371b666dd 100644
--- a/OpenPNM/Geometry/models/throat_offset_vertices.py
+++ b/OpenPNM/Geometry/models/throat_offset_vertices.py
@@ -5,40 +5,9 @@ throat_offset_vertices -- Offeset throat vertices using a fibre radius parameter
 
 """
 import scipy as sp
-import OpenPNM.Utilities.vertexops as vo
 import OpenPNM.Utilities.transformations as tr
-
-
-def voronoi(network, geometry, offset, **kwargs):
-    r"""
-    Offset the throat vertices effectively erroding the facet by the offset distance
-    supplied.
-    """
-    Nt = geometry.num_throats()
-    area = sp.ndarray(Nt)
-    perimeter = sp.ndarray(Nt)
-    offset_verts = sp.ndarray(Nt, dtype=object)
-    offset_error = sp.ndarray(Nt)
-    throat_COM = sp.ndarray([Nt, 3])
-    for i in range(Nt):
-        offset_rand = (sp.random.rand(1)-0.5)*offset
-        throat_verts = geometry['throat.vertices'][i]
-        throat_normal = geometry['throat.normal'][i]
-        area[i], perimeter[i], offset_verts[i], throat_COM[i], offset_error[i] = \
-            vo.get_throat_geom(throat_verts, throat_normal, offset)
-
-    for i in range(Nt):
-        if offset_error[i] > 0 and len(offset_verts[i]) > 0:
-            offset_verts[i] = []
-    # Properties that depend on the offset vertices are the area, perimeter and the
-    # centroid or COM. To speed things up we could save them all now rather than
-    # processing them individually.
-    if kwargs['set_dependent'] is True:
-        geometry['throat.area'] = area
-        geometry['throat.perimeter'] = perimeter
-        geometry['throat.centroid'] = throat_COM
-
-    return offset_verts
+from OpenPNM.Base import logging
+logger = logging.getLogger(__name__)
 
 
 def distance_transform(network, geometry, offset, **kwargs):
@@ -67,6 +36,7 @@ def distance_transform(network, geometry, offset, **kwargs):
     z_axis = [0, 0, 1]
 
     for i in range(Nt):
+        logger.info("Processing throat " + str(i+1)+" of "+str(Nt))
         # For boundaries some facets will already be aligned with the axis - if this
         # is the case a rotation is unnecessary and could also cause problems
         angle = tr.angle_between_vectors(normals[i], z_axis)
@@ -115,6 +85,9 @@ def distance_transform(network, geometry, offset, **kwargs):
             z_plane = sp.unique(np.around(z, order+2))
             if len(z_plane) > 1:
                 print('Rotation for image analysis failed')
+                temp_arr = np.ones(1)
+                temp_arr.fill(np.mean(z_plane))
+                z_plane = temp_arr
             "Fill in the convex hull polygon"
             convhullimg = convex_hull_image(img_pad)
             # Perform a Distance Transform and black out points less than r to create
diff --git a/OpenPNM/Geometry/models/throat_perimeter.py b/OpenPNM/Geometry/models/throat_perimeter.py
index d815b18b0..dbfbe57d1 100644
--- a/OpenPNM/Geometry/models/throat_perimeter.py
+++ b/OpenPNM/Geometry/models/throat_perimeter.py
@@ -7,6 +7,7 @@ Submodule -- throat_perimeter
 import scipy as _sp
 import OpenPNM.Utilities.transformations as tr
 import OpenPNM.Utilities.vertexops as vo
+from scipy.spatial import ConvexHull
 
 
 def voronoi(geometry, **kwargs):
@@ -20,6 +21,9 @@ def voronoi(geometry, **kwargs):
     for i in range(Nt):
         if len(verts[i]) > 2:
             verts_2D = tr.rotate_and_chop(verts[i], normals[i], [0, 0, 1])
+            # Get in hull order
+            hull = ConvexHull(verts_2D, qhull_options='QJ Pp')
+            verts_2D = verts_2D[hull.vertices]
             perimeter[i] = vo.PolyPerimeter2D(verts_2D)
         else:
             perimeter[i] = 0.0
diff --git a/OpenPNM/Geometry/models/throat_shape_factor.py b/OpenPNM/Geometry/models/throat_shape_factor.py
index 6aa924aca..79bbf6a33 100644
--- a/OpenPNM/Geometry/models/throat_shape_factor.py
+++ b/OpenPNM/Geometry/models/throat_shape_factor.py
@@ -23,15 +23,20 @@ def compactness(geometry, throat_perimeter='throat.perimeter',
     C = _sp.ones(geometry.num_throats())
     C[ts] = P[ts]**2/A[ts]
     verts = geometry['throat.offset_vertices']
-    for i in range(len(verts)):
+    alpha = _sp.ones_like(C)
+    for i in ts:
         if len(verts[i]) == 3:
             # Triangular Correction
-            C[i] = C[i]*(25/17) + (40*_sp.sqrt(3)/17)
+            alpha[i] = C[i]*(25/17) + (40*_sp.sqrt(3)/17)
         elif len(verts[i]) == 4:
             # Rectangular Correction
-            C[i] = C[i]*(22/7) + (65/3)
+            alpha[i] = C[i]*(22/7) - (65/3)
         elif len(verts[i]) > 4:
             # Approximate Elliptical Correction
-            C[i] = C[i]*(8/3) + (8*_sp.pi/3)
+            alpha[i] = C[i]*(8/3) - (8*_sp.pi/3)
+    # For a perfect circle alpha = 8*pi so normalize by this
+    alpha /= 8*_sp.pi
+    # Very small throats could have values less than one
+    alpha[alpha < 1.0] = 1.0
 
-    return C
+    return alpha
diff --git a/OpenPNM/Network/__Cubic__.py b/OpenPNM/Network/__Cubic__.py
index 37bb0248a..d194fd984 100644
--- a/OpenPNM/Network/__Cubic__.py
+++ b/OpenPNM/Network/__Cubic__.py
@@ -43,7 +43,7 @@ class Cubic(GenericNetwork):
         - 26: Faces, Edges and Corners
 
     template : array of booleans
-        An (i,j,k) array with True where the Network should be defiend and
+        An (i,j,k) array with True where the Network should be defined and
         False elsewhere. This approach is useful for creating networks of non-
         cuboid shape like spheres or cylinders, but still with a cubic lattice
         topology.
diff --git a/OpenPNM/Network/__DelaunayCubic__.py b/OpenPNM/Network/__DelaunayCubic__.py
index 604bd2f70..ac7683146 100644
--- a/OpenPNM/Network/__DelaunayCubic__.py
+++ b/OpenPNM/Network/__DelaunayCubic__.py
@@ -27,6 +27,19 @@ class DelaunayCubic(Delaunay):
     name : string
         A unique name for the network
 
+    shape : tuple of ints
+        The (i,j,k) size and shape of the network.
+
+    spacing : 3 x 1 array defining the base lattice spacing of the network
+
+    perturbation : float between 0 and 1 controlling the maximum perturbation
+        of lattice points as a fraction of the lattice spacing
+
+    arrangement : string
+        usage:  'SC'- Simple Cubic (default if left blank)
+                'O' - Orthorhombic
+                'BCC' - Body Centred Cubic
+                'FCC' - Face Centred Cubic
     Examples
     --------
     >>> import OpenPNM
@@ -38,12 +51,10 @@ class DelaunayCubic(Delaunay):
 
     """
 
-    def __init__(self, shape=None, template=None, spacing=[1, 1, 1],
-                 jiggle_factor=0.1, arrangement='SC', **kwargs):
+    def __init__(self, shape=None, spacing=[1, 1, 1],
+                 perturbation=0.1, arrangement='SC', **kwargs):
         if shape is not None:
             self._arr = np.atleast_3d(np.empty(shape))
-        elif template is not None:
-            self._arr = sp.array(template, ndmin=3, dtype=bool)
         else:
             self._arr = np.atleast_3d(np.empty([3, 3, 3]))
 
@@ -53,13 +64,13 @@ class DelaunayCubic(Delaunay):
         self._spacing = sp.asarray(spacing)
         self._num_pores = np.prod(np.asarray(self._shape))
         self._domain_size = np.asarray(self._shape) * self._spacing
-        self._jiggle_factor = jiggle_factor
+        self._perturbation = perturbation
         self._arrangement = arrangement
         super().__init__(num_pores=self._num_pores,
                          domain_size=self._domain_size,
                          **kwargs)
 
-    def _generate_pores(self):
+    def _generate_pores(self, prob=None):
         r"""
         Generate the pores with numbering scheme.
         """
@@ -116,7 +127,7 @@ class DelaunayCubic(Delaunay):
             face_points = np.asarray(face_points)
             points = np.concatenate((points, face_points))
 
-        jiggle = (np.random.rand(len(points), 3)-0.5)*self._jiggle_factor
+        jiggle = (np.random.rand(len(points), 3)-0.5)*self._perturbation
         points += jiggle
         points *= self._spacing
 
diff --git a/OpenPNM/Network/__Delaunay__.py b/OpenPNM/Network/__Delaunay__.py
index ac0d4e205..1580f1f5b 100644
--- a/OpenPNM/Network/__Delaunay__.py
+++ b/OpenPNM/Network/__Delaunay__.py
@@ -13,8 +13,6 @@ import scipy.sparse as sprs
 import scipy.spatial as sptl
 import scipy.ndimage as spim
 from scipy.spatial import Voronoi
-from scipy import stats as st
-from scipy.special import cbrt
 from OpenPNM.Network import GenericNetwork
 from OpenPNM.Base import logging
 logger = logging.getLogger(__name__)
@@ -25,13 +23,20 @@ class Delaunay(GenericNetwork):
     This class contains the methods for creating a *Delaunay* network topology
     based connecting pores with a Delaunay tessellation.
 
-    To invoke the actual generation it is necessary to run the `generate` method.
-
     Parameters
     ----------
     name : string
         A unique name for the network
-
+    domain_size : list of floats, [Lx,Ly,Lz]
+        Bounding cube for internal pore positions
+    num_pores : int
+        Number of pores to place randomly within domain
+    prob : 3D float array
+        Values should be between 0 and 1 as determines probability of point
+        with relative domain coordinates being kept. Array does not have to
+        be same size as domain because positions are re-scaled
+    base_points : [Np,3] float array
+        coordinates to use instead of random generation
     Examples
     --------
     >>> import OpenPNM
@@ -42,91 +47,80 @@ class Delaunay(GenericNetwork):
 
     """
 
-    def __init__(self, num_pores=None, domain_size=None, **kwargs):
+    def __init__(self, num_pores=None, domain_size=None, prob=None,
+                 base_points=None, **kwargs):
         """
         Create Delauny network object
         """
         super().__init__(**kwargs)
-        if (num_pores and domain_size) is None:
-            num_pores = 1
-            domain_size = [1.0, 1.0, 1.0]
-        else:
-            self.generate(num_pores, domain_size)
+        self.generate(num_pores, domain_size, prob, base_points)
 
-    def generate(self, num_pores, domain_size):
+    def generate(self, num_pores, domain_size, prob, base_points):
         r"""
         Method to trigger the generation of the network
-
-        Parameters
-        ----------
-        domain_size : list of floats, [Lx,Ly,Lz]
-            Bounding cube for internal pore positions
-        num_pores : int
-            Number of pores to place randomly within domain
-
         """
         logger.info('Start of network topology generation')
-        self._generate_setup(num_pores, domain_size)
-        self._generate_pores()
+        self._generate_setup(num_pores, domain_size, base_points)
+        if base_points is not None:
+            try:
+                dim = sp.shape(base_points)[1]
+                if dim != 3:
+                    raise Exception('base points must be 3D')
+            except:
+                raise Exception('base points must be 3D')
+            self['pore.coords'] = base_points
+        else:
+            self._generate_pores(prob)
         self._generate_throats()
         logger.debug('Network generation complete')
 
-    def _generate_setup(self, num_pores, domain_size):
+    def _generate_setup(self, num_pores, domain_size, base_points):
         r"""
         Perform applicable preliminary checks and calculations required for
         generation
         """
-        logger.debug('generate_setup: Perform preliminary calculations')
-        if domain_size is not None and num_pores is not None:
-            self._Lx = domain_size[0]
-            self._Ly = domain_size[1]
-            self._Lz = domain_size[2]
-            self._Np = num_pores
-            r"""
-            TODO: Fix this, btype should be received as an argument
-            """
-            self._btype = [0, 0, 0]
-        else:
-            logger.error('domain_size and num_pores must be specified')
-            raise Exception('domain_size and num_pores must be specified')
+        logger.debug('generate_setup: Perform preliminary calculations and checks')
+        if domain_size is None:
+            raise Exception('domain_size must always be specified')
+        if num_pores is None and base_points is None:
+            raise Exception('num_pores or base_points must be specified')
+        elif num_pores is None and base_points is not None:
+            num_pores = len(base_points)
+        elif num_pores is not None and base_points is not None:
+            logger.warning('both num_pores and base_points arguments given' +
+                           ' num_pores over-written')
+            num_pores = len(base_points)
+
+        self._Lx = domain_size[0]
+        self._Ly = domain_size[1]
+        self._Lz = domain_size[2]
+        self._Np = num_pores
+        r"""
+        TODO: Fix this, btype should be received as an argument
+        """
+        self._btype = [0, 0, 0]
 
-    def _generate_pores(self):
+    def _generate_pores(self, prob=None):
         r"""
         Generate the pores with numbering scheme.
         """
         logger.info('Place randomly located pores in the domain')
-        coords = np.zeros([self._Np, 3])
-        # Reject points close to boundaries - if False there will be slightly more
-        rejection = [False, False, True]
-        for j in range(3):
+        if prob is not None:
+            coords = []
             i = 0
             while i < self._Np:
-                coord = np.random.uniform(0, 1, 1)
-                if self._reject(coord) == rejection[j]:
-                    coords[i][j] = coord
+                coord = np.random.rand(3)
+                [indx, indy, indz] = np.floor(coord*np.shape(prob)).astype(int)
+                p = prob[indx][indy][indz]
+                if np.random.rand(1) <= p:
+                    coords.append(coord)
                     i += 1
-        coords *= np.array([self._Lx, self._Ly, self._Lz])
-
-        self['pore.coords'] = coords
-
-    def _prob_func(self, m):
-        a = 35
-        b = 0.2
-        p = ((m**a) + ((1-m)**a) + (2*b))/(1 + (2*b))
-        return p
-
-    def _reject(self, point):
-
-        P = self._prob_func(point)
-        nrand = np.random.uniform(0, 1, 1)
-        # Place more points at the sides of the domain and fewer at the
-        # top and bottom
-        if P < nrand:
-            rejection = True
+            coords = np.asarray(coords)
         else:
-            rejection = False
+            coords = np.random.random([self._Np, 3])
 
-        return rejection
+        coords *= np.array([self._Lx, self._Ly, self._Lz])
+        self['pore.coords'] = coords
 
     def _generate_throats(self):
         r"""
@@ -218,248 +212,6 @@ class Delaunay(GenericNetwork):
         self['throat.vert_index'] = throat_verts
         logger.debug(sys._getframe().f_code.co_name + ': End of method')
 
-    def _add_labels(self):
-        r"""
-        Deprecated if using add_boundaries()
-        This finds surface pores simply by proximity to the domain boundaries.
-        A better approach is necessary
-        """
-        coords = self['pore.coords']
-        self['pore.front'] = coords[:, 0] < (0.1*self._Lx)
-        self['pore.back'] = coords[:, 0] > (0.9*self._Lx)
-        self['pore.left'] = coords[:, 1] < (0.1*self._Ly)
-        self['pore.right'] = coords[:, 1] > (0.9*self._Ly)
-        self['pore.bottom'] = coords[:, 2] < (0.1*self._Lz)
-        self['pore.top'] = coords[:, 2] > (0.9*self._Lz)
-        bnds = self.pores(labels=['front', 'back', 'left', 'right', 'bottom', 'top'])
-        self['pore.boundary'] = False
-        self['pore.boundary'] = bnds
-
-    def _add_boundaries(self):
-        r"""
-        This is an alternative means of adding boundaries
-        """
-        logger.info('add_boundaries: start of method')
-
-        import scipy.spatial as sptl
-        import scipy.sparse as sprs
-        Lx = self._Lx
-        Ly = self._Ly
-        Lz = self._Lz
-        Np = self.num_pores()
-        btype = self._btype
-        boffset = 0.05
-
-        # Translate internal pores to each face of domain
-        poffset = np.zeros((7, 3))
-        poffset[[2, 5], 0] = [-Lx, Lx]
-        poffset[[3, 4], 1] = [-Ly, Ly]
-        poffset[[1, 6], 2] = [-Lz, Lz]
-        pcoords = pcoords0 = self['pore.coords']
-        for i in np.r_[1:7]:
-            pcoords = np.concatenate((pcoords, pcoords0 + poffset[i, :]), axis=0)
-
-        # Use some twisted logic to get bval list of + for boundary and -
-        # for periodic faces
-        bval = [0, 1, 2, 3, 4, 5, 6] * \
-            (np.array([0, btype[2], btype[0], btype[1],
-                       btype[1], btype[0], btype[2]])*-2+1)
-        ptype = np.zeros((Np,), dtype=int)
-        for i in np.r_[1:7]:
-            ptype = \
-                np.concatenate((ptype, np.ones((Np,), dtype=int)*bval[i]), axis=0)
-
-        # pnum contains the internal ID number of the boundary pores
-        # for connecting periodic points
-        pnum = self.pores()
-        pnum = np.tile(pnum, 7)
-
-        Tri = sptl.Delaunay(pcoords)
-        adjmat = \
-            sprs.lil_matrix((np.shape(pcoords)[0], np.shape(pcoords)[0]), dtype=int)
-        for i in np.arange(0, np.shape(Tri.simplices)[0]):
-            # Keep only simplices that are fully in real domain
-            adjmat[Tri.simplices[i], Tri.simplices[i]] = 1
-        adjmat = sprs.triu(adjmat, k=1, format='lil')
-        for i in np.arange(0, Np):
-            # Add periodic throats to the netowrk (if any)
-            tpore2 = pnum[adjmat.rows[i]][ptype[adjmat.rows[i]] < 0]
-            tpore1 = np.ones_like(tpore2, dtype=int) * i
-            conns = self['throat.conns']
-            conns = np.concatenate((conns, np.vstack((tpore1, tpore2)).T), axis=0)
-            # Add boundary pores and throats to the network
-            newporetyps = np.unique(ptype[adjmat.rows[i]][ptype[adjmat.rows[i]] > 0])
-            newporenums = \
-                np.r_[self.num_pores():self.num_pores()+np.size(newporetyps)]
-            tpore2 = newporenums
-            tpore1 = np.ones_like(tpore2, dtype=int) * i
-            conns = np.concatenate((conns, np.vstack((tpore1, tpore2)).T), axis=0)
-            self['throat.conns'] = conns
-            bcoords = np.zeros((7, 3), dtype=float)
-            coords = self['pore.coords']
-            bcoords[1, :] = [coords[i, 0], coords[i, 1], 0-Lz*boffset]
-            bcoords[2, :] = [0-Lx*boffset, coords[i, 1], coords[i, 2]]
-            bcoords[3, :] = [coords[i, 0], -Ly*boffset, coords[i, 2]]
-            bcoords[4, :] = [coords[i, 0], Ly+Ly*boffset, coords[i, 2]]
-            bcoords[5, :] = [Lx+Lx*boffset, coords[i, 1], coords[i, 2]]
-            bcoords[6, :] = [coords[i, 0], coords[i, 1], Lz+Lz*boffset]
-            newporecoords = bcoords[newporetyps, :]
-            coords = np.concatenate((coords, newporecoords), axis=0)
-            self['pore.coords'] = coords
-        # Reset number of pores and throats (easier than tracking it)
-        nums = np.r_[0:np.shape(coords)[0]]
-        self['pore.numbering'] = nums
-        self['pore.numbering'] = np.ones((nums[-1]+1,), dtype=bool)
-        nums = np.r_[0:np.shape(conns)[0]]
-        self['throat.numbering'] = nums
-        self['throat.numbering'] = np.ones((nums[-1]+1,), dtype=bool)
-        logger.debug('add_boundaries: end of method')
-
-    def _add_boundaries_old(self):
-        logger.info('add_boundaries_old: Start of method')
-
-        self.add_opposing_boundaries(btype=[2, 5])
-        self.add_opposing_boundaries(btype=[3, 4])
-        self.add_opposing_boundaries(btype=[1, 6])
-
-    def _add_opposing_boundaries(self, btype=[1, 6]):
-        r"""
-        btype indicates which two boundaries are being added by type
-        """
-        logger.info('add_opposing_boundaries: start of method')
-
-        if btype == [2, 5]:
-            D = 0
-            W = 1
-            H = 2
-        elif btype == [3, 4]:
-            D = 1
-            W = 0
-            H = 2
-        elif btype == [1, 6]:
-            D = 2
-            W = 1
-            H = 0
-
-        Lx = self.domain_size[D]
-        Ly = self.domain_size[W]
-        Lz = self.domain_size[H]
-        # Rotate pore coordinates (use only internal pores)
-        pnum = self._net.pore_data['numbering'][self._net.pore_data['type'] == 0]
-        pcoords = np.zeros_like(self._net.pore_data['coords'][pnum, :])
-        pcoords[:, 0] = self._net.pore_data['coords'][pnum, D]
-        pcoords[:, 1] = self._net.pore_data['coords'][pnum, W]
-        pcoords[:, 2] = self._net.pore_data['coords'][pnum, H]
-
-        # Determine dimensions of image from dimensions of domain
-        f = 100  # minimum image dimension
-        im_dim = [0, 0, 0]
-        im_dim[0] = np.floor(f*Lx/np.min([Lx, Ly, Lz]))
-        im_dim[1] = np.floor(f*Ly/np.min([Lx, Ly, Lz]))
-        im_dim[2] = np.floor(f*Lz/np.min([Lx, Ly, Lz]))
-        im_dim = np.array(im_dim, dtype=int)
-
-        # Convert pore coordinates into image subscripts
-        im_subs = np.zeros_like(pcoords, dtype=int)
-        im_subs[:, 0] = pcoords[:, 0]*im_dim[0]/Lx
-        im_subs[:, 1] = pcoords[:, 1]*im_dim[1]/Ly
-        im_subs[:, 2] = pcoords[:, 2]*im_dim[2]/Lz
-        # Find linear indices of each pore in the new image
-        im_inds = np.ravel_multi_index((im_subs[:, 0], im_subs[:, 1], im_subs[:, 2]),
-                                       dims=(im_dim), order='F')
-
-        # Generate 3D image of points (place pore numbers at each site for use later)
-        img = np.zeros(im_dim, dtype=int)
-        img[im_subs[:, 0], im_subs[:, 1], im_subs[:, 2]] = pnum
-
-        # Perform distance transform on points and also get 'indicies' of each point
-        img_dt, ind_dt = spim.distance_transform_edt(img == 0)
-
-        # Project all* internal points to x face
-        # *Note that it's possible/likely that mutliple internal points map to the
-        # same boundary point
-        img_bd0 = np.zeros([im_dim[1], im_dim[2]], dtype=int)
-        img_bd1 = np.zeros([im_dim[1], im_dim[2]], dtype=int)
-        img_bd0[im_subs[:, 1], im_subs[:, 2]] = im_inds
-        img_bd1[im_subs[:, 1], im_subs[:, 2]] = im_inds
-
-        # Create 2D array of distance transform indices for 0 and end faces
-        dt_D0 = ind_dt[0, 0, :, :]*(img_bd0 > 0)  # 0 face
-        dt_D1 = ind_dt[0, -1, :, :]*(img_bd1 > 0)  # end face
-
-        # Create a 2D mask containing x coordinates of internal points
-        # and -1 elsewhere
-        img_D0 = -np.ones([im_dim[1], im_dim[2]], dtype=int)
-        img_D1 = -np.ones([im_dim[1], im_dim[2]], dtype=int)
-        img_D0[im_subs[:, 1], im_subs[:, 2]] = im_subs[:, 0]
-        img_D1[im_subs[:, 1], im_subs[:, 2]] = im_subs[:, 0]
-
-        # Find where x value of internal points corresponds to x value of
-        # distance transform indices
-        img_bd0 = (img_D0 == dt_D0)*img_bd0
-        img_bd1 = (img_D1 == dt_D1)*img_bd1
-
-        # Convert boundary sites to linear indices
-        inds_bd0 = img_bd0[np.nonzero(img_bd0)]
-        inds_bd1 = img_bd1[np.nonzero(img_bd1)]
-
-        # Use linear indices to find pore ID nums
-        nums_bd0 = img[np.unravel_index(inds_bd0, dims=(im_dim), order='F')]
-        nums_bd1 = img[np.unravel_index(inds_bd1, dims=(im_dim), order='F')]
-        nums_bd = np.append(nums_bd0, nums_bd1)
-        types_bd = np.append(np.zeros_like(nums_bd0), np.ones_like(nums_bd1))
-
-        # Add new boundary pores and throats to the network
-        # Get all pores including previously added boundaries
-        Np = self._net.num_pores()
-        bp_numbering = np.r_[Np:Np+np.size(nums_bd)]
-        bp_type = (types_bd == 0)*btype[0] + (types_bd == 1)*btype[1]
-        bp_coords = np.zeros([np.size(nums_bd), 3])
-        bp_coords[types_bd == 0, D] = np.zeros_like(nums_bd0) - .0001
-        bp_coords[types_bd == 0, W] = pcoords[nums_bd0, 1]
-        bp_coords[types_bd == 0, H] = pcoords[nums_bd0, 2]
-        bp_coords[types_bd == 1, D] = np.ones_like(nums_bd1)*Lx+0.0001
-        bp_coords[types_bd == 1, W] = pcoords[nums_bd1, 1]
-        bp_coords[types_bd == 1, H] = pcoords[nums_bd1, 2]
-        self._net.pore_data['numbering'] = \
-            np.append(self._net.pore_data['numbering'], bp_numbering)
-        self._net.pore_data['type'] = np.append(self._net.pore_data['type'], bp_type)
-        self._net.pore_data['coords'] = \
-            np.concatenate((self._net.pore_data['coords'], bp_coords))
-        Nt = self._net.num_throats()
-        bt_numbering = np.r_[Nt:Nt + np.size(nums_bd)]
-        bt_type = np.ones(np.size(nums_bd), dtype=int)*2
-        bt_connections = np.zeros([np.size(nums_bd), 2], dtype=int)
-        bt_connections[:, 0] = nums_bd
-        bt_connections[:, 1] = bp_numbering
-        self._net.throat_data['numbering'] = \
-            np.append(self._net.throat_data['numbering'], bt_numbering)
-        self._net.throat_data['type'] = \
-            np.append(self._net.throat_data['type'], bt_type)
-        self._net.throat_data['conns'] = \
-            np.concatenate((self._net.throat_data['conns'], bt_connections))
-
-    def domain_size(self, dimension=''):
-        r"""
-        This is a simple way to find the domain sizes.
-        N.B
-        Will not work with saved and loaded networks
-        """
-        if dimension == 'front' or dimension == 'back':
-            return self._Ly*self._Lz
-        if dimension == 'left' or dimension == 'right':
-            return self._Lx*self._Lz
-        if dimension == 'top' or dimension == 'bottom':
-            return self._Lx*self._Ly
-        if dimension == 'volume':
-            return self._Lx*self._Ly*self._Lz
-        if dimension == 'height':
-            return self._Lz
-        if dimension == 'width':
-            return self._Lx
-        if dimension == 'depth':
-            return self._Ly
-
     def add_boundaries(self):
 
         r"""
@@ -535,7 +287,8 @@ class Delaunay(GenericNetwork):
                     elif len(np.unique(throat_verts[:, 2])) == 1:
                         new_pore_coord[2] = np.unique(throat_verts[:, 2])
                     else:
-                        new_pore_coord = throat_verts.mean()
+                        new_pore_coord = np.mean(throat_verts, axis=0)
+                        pass
                     bound_coords.append(new_pore_coord)
                     bound_conns.append(np.array([my_pore, new_throat_count + Np]))
                     bound_vert_index.append(dict(zip(v, throat_verts)))
@@ -556,9 +309,14 @@ class Delaunay(GenericNetwork):
         right = self.pores()[self['pore.coords'][:, 1] == max_point[1]]
         bottom = self.pores()[self['pore.coords'][:, 2] == min_point[2]]
         top = self.pores()[self['pore.coords'][:, 2] == max_point[2]]
+        if len(top) == 0:
+            top = self.pores()[self['pore.coords'][:, 2] ==
+                               np.asarray(bound_coords)[:, 2].max()]
         # Assign labels
         self['pore.boundary'] = False
         self['pore.boundary'][new_pore_ids] = True
+        self['throat.boundary'] = False
+        self['throat.boundary'][new_throat_ids] = True
         self['pore.right_boundary'] = False
         self['pore.left_boundary'] = False
         self['pore.front_boundary'] = False
@@ -579,16 +337,6 @@ class Delaunay(GenericNetwork):
         r"""
         Returns the distance between two faces
         No coplanar checking this is done in vertex_dimension
-
-        Example
-        --------
-        >>> import OpenPNM
-        >>> pn = OpenPNM.Network.Delaunay(num_pores=100, domain_size=[3,2,1])
-        >>> pn.add_boundaries()
-        >>> B1 = pn.pores('left_boundary')
-        >>> B2 = pn.pores('right_boundary')
-        >>> pn.domain_length(B1,B2)
-        2.0
         """
         L = vo.vertex_dimension(self, face_1, face_2, parm='length')
         return L
@@ -597,15 +345,6 @@ class Delaunay(GenericNetwork):
         r"""
         Returns the area of a face
         No coplanar checking this is done in vertex_dimension
-        Example
-        --------
-        >>> import OpenPNM
-        >>> pn = OpenPNM.Network.Delaunay(num_pores=100, domain_size=[3,2,1])
-        >>> pn.add_boundaries()
-        >>> B1 = pn.pores('left_boundary')
-        >>> B2 = pn.pores('right_boundary')
-        >>> pn.domain_area(B1)
-        3.0
         """
         A = vo.vertex_dimension(self, face, parm='area')
 
@@ -616,10 +355,41 @@ class Delaunay(GenericNetwork):
         After the offsetting routine throats with zero area have been fully occluded.
         Remove these from the network and also remove pores that are isolated
         """
-        occluded_ts = list(self.throats()[self['throat.area'] == 0])
-        if len(occluded_ts) > 0:
+        occluded_ts = self['throat.area'] == 0
+        if np.sum(occluded_ts) > 0:
             self.trim(throats=occluded_ts)
         # Also get rid of isolated pores
         isolated_ps = self.check_network_health()['isolated_pores']
         if len(isolated_ps) > 0:
             self.trim(isolated_ps)
+
+    def export_vor_fibres(self):
+        r"""
+        Run through the throat vertices, compute the convex hull order and save
+        the vertices and ordered faces in a pickle dictionary to be used in
+        blender
+        """
+        import pickle as pickle
+        Indices = []
+        for t in self.throats():
+            indices = list(self["throat.vert_index"][t].keys())
+            verts = self._vor.vertices[indices]
+            # Need to order the indices in convex hull order
+            # Compute the standard deviation in all coordinates and eliminate
+            # the axis with the smallest to make 2d
+            stds = [np.std(verts[:, 0]), np.std(verts[:, 1]), np.std(verts[:, 2])]
+            if np.argmin(stds) == 0:
+                verts2d = np.vstack((verts[:, 1], verts[:, 2])).T
+            elif np.argmin(stds) == 1:
+                verts2d = np.vstack((verts[:, 0], verts[:, 2])).T
+            else:
+                verts2d = np.vstack((verts[:, 0], verts[:, 1])).T
+            # 2d convexhull returns vertices in hull order
+            hull2d = sptl.ConvexHull(verts2d, qhull_options='QJ Pp')
+            # Re-order the vertices and save as list (blender likes them as lists)
+            Indices.append(np.asarray(indices)[hull2d.vertices].tolist())
+        # Create dictionary to pickle
+        data = {}
+        data["Verts"] = self._vor.vertices
+        data["Indices"] = Indices
+        pickle.dump(data, open("fibres.p", "wb"))
diff --git a/OpenPNM/Physics/__GenericPhysics__.py b/OpenPNM/Physics/__GenericPhysics__.py
index bc5209117..fb7e2d1ac 100644
--- a/OpenPNM/Physics/__GenericPhysics__.py
+++ b/OpenPNM/Physics/__GenericPhysics__.py
@@ -98,10 +98,10 @@ class GenericPhysics(OpenPNM.Base.Core):
             Indicates whether list of pores or throats is to be added or removed
             from the object.  Options are 'add' (default) or 'remove'.
         """
-        if len(pores) > 0:
+        if sp.size(pores) > 0:
             pores = sp.array(pores, ndmin=1)
             self._set_locations(element='pore', locations=pores, mode=mode)
-        if len(throats) > 0:
+        if sp.size(throats) > 0:
             throats = sp.array(throats, ndmin=1)
             self._set_locations(element='throat', locations=throats, mode=mode)
 
diff --git a/OpenPNM/Physics/models/diffusive_conductance.py b/OpenPNM/Physics/models/diffusive_conductance.py
index b88fab87e..2544e436a 100644
--- a/OpenPNM/Physics/models/diffusive_conductance.py
+++ b/OpenPNM/Physics/models/diffusive_conductance.py
@@ -13,7 +13,7 @@ def bulk_diffusion(physics, phase, network, pore_molar_density='pore.molar_densi
                    pore_diffusivity='pore.diffusivity', pore_area='pore.area',
                    pore_diameter='pore.diameter', throat_area='throat.area',
                    throat_length='throat.length', throat_diameter='throat.diameter',
-                   calc_pore_len=True, **kwargs):
+                   calc_pore_len=True, shape_factor='throat.shape_factor', **kwargs):
     r"""
     Calculate the diffusive conductance of conduits in network, where a
     conduit is ( 1/2 pore - full throat - 1/2 pore ) based on the areas
@@ -40,7 +40,8 @@ def bulk_diffusion(physics, phase, network, pore_molar_density='pore.molar_densi
     parea = network[pore_area]
     pdia = network[pore_diameter]
     # Get the properties of every throat
-    tarea = network[throat_area]
+    tdia = network[throat_diameter]
+    tarea = _sp.pi*(tdia/2)**2
     tlen = network[throat_length]
     # Interpolate pore phase property values to throats
     cp = phase[pore_molar_density]
@@ -67,7 +68,14 @@ def bulk_diffusion(physics, phase, network, pore_molar_density='pore.molar_densi
     gp2[~(gp2 > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
     # Find g for full throat, remove any non-positive lengths
     tlen[tlen <= 0] = 1e-12
-    gt = ct*DABt*tarea/tlen
+    # Get shape factor
+    try:
+        sf = network[shape_factor]
+    except:
+        sf = _sp.ones(network.num_throats())
+    gt = (1/sf)*ct*DABt*tarea/tlen
+    # Set 0 conductance pores (boundaries) to inf
+    gt[~(gt > 0)] = _sp.inf
     value = (1/gt + 1/gp1 + 1/gp2)**(-1)
     value = value[phase.throats(physics.name)]
     return value
diff --git a/OpenPNM/Physics/models/hydraulic_conductance.py b/OpenPNM/Physics/models/hydraulic_conductance.py
index 9f23fb7e0..4ed5c0a47 100644
--- a/OpenPNM/Physics/models/hydraulic_conductance.py
+++ b/OpenPNM/Physics/models/hydraulic_conductance.py
@@ -12,6 +12,7 @@ import OpenPNM.Utilities.misc as misc
 def hagen_poiseuille(physics, phase, network, pore_diameter='pore.diameter',
                      pore_viscosity='pore.viscosity', throat_length='throat.length',
                      throat_diameter='throat.diameter', calc_pore_len=True,
+                     shape_factor='throat.shape_factor',
                      **kwargs):
     r"""
     Calculates the hydraulic conductivity of throat assuming cylindrical
@@ -62,7 +63,12 @@ def hagen_poiseuille(physics, phase, network, pore_diameter='pore.diameter',
     tlen = network[throat_length]
     # Remove any non-positive lengths
     tlen[tlen <= 0] = 1e-12
-    gt = _sp.pi*(tdia)**4/(128*tlen*mut)
+    # Get shape factor
+    try:
+        sf = network[shape_factor]
+    except:
+        sf = _sp.ones(network.num_throats())
+    gt = (1/sf)*_sp.pi*(tdia)**4/(128*tlen*mut)
     gt[~(gt > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
     value = (1/gt + 1/gp1 + 1/gp2)**(-1)
     value = value[phase.throats(physics.name)]
diff --git a/OpenPNM/Physics/models/multiphase.py b/OpenPNM/Physics/models/multiphase.py
index 69e1e5693..d11a20b81 100644
--- a/OpenPNM/Physics/models/multiphase.py
+++ b/OpenPNM/Physics/models/multiphase.py
@@ -50,20 +50,20 @@ def conduit_conductance(physics, phase, network, throat_conductance,
     """
     throats = phase.throats(physics.name)
     if mode == 'loose':
-        closed_conduits = -sp.array(phase[throat_occupancy], dtype=bool)
+        closed_conduits = ~sp.array(phase[throat_occupancy], dtype=bool)
     else:
-        throats_closed = -sp.array(phase[throat_occupancy], dtype=bool)
+        throats_closed = ~sp.array(phase[throat_occupancy], dtype=bool)
         connected_pores = network.find_connected_pores(throats)
         pores_1 = connected_pores[:, 0]
         pores_2 = connected_pores[:, 1]
-        pores_1_closed = -sp.array(phase[pore_occupancy][pores_1], dtype=bool)
-        pores_2_closed = -sp.array(phase[pore_occupancy][pores_2], dtype=bool)
+        pores_1_closed = ~sp.array(phase[pore_occupancy][pores_1], dtype=bool)
+        pores_2_closed = ~sp.array(phase[pore_occupancy][pores_2], dtype=bool)
         if mode == 'medium':
             closed_conduits = throats_closed | (pores_1_closed & pores_2_closed)
 
         if mode == 'strict':
             closed_conduits = pores_1_closed | throats_closed | pores_2_closed
-    open_conduits = -closed_conduits
+    open_conduits = ~closed_conduits
     throat_value = phase[throat_conductance]
     value = throat_value*open_conduits + throat_value*closed_conduits*factor
     return value
@@ -99,7 +99,10 @@ def late_pore_filling(physics, phase, network, Pc, Swp_star=0.2, eta=3,
     prop = phase[throat_capillary_pressure]
     neighborTs = network.find_neighbor_throats(pores, flatten=False)
     Pc_star = sp.array([sp.amin(prop[row]) for row in neighborTs])
-    Swp = Swp_star*(Pc_star/Pc)**eta
+    if Pc > 0:
+        Swp = Swp_star*(Pc_star/Pc)**eta
+    else:
+        Swp = sp.zeros(len(Pc_star))
     if wetting_phase:
         values = Swp*phase[pore_occupancy]*(Pc_star < Pc)
     else:
diff --git a/OpenPNM/Postprocessing/Plots.py b/OpenPNM/Postprocessing/Plots.py
index eab2be8a0..a2a35f2ab 100644
--- a/OpenPNM/Postprocessing/Plots.py
+++ b/OpenPNM/Postprocessing/Plots.py
@@ -209,6 +209,11 @@ def distributions(net, fig=None, throat_diameter='throat.diameter',
     else:
         include_pores = net['pore.all']
         include_throats = net['throat.all']
+
+    if exclude_boundaries:
+        include_pores *= ~net["pore.boundary"]
+        include_throats *= ~net["throat.boundary"]
+
     pores = net.pores()[include_pores]
     throats = net.throats()[include_throats]
 
diff --git a/OpenPNM/Utilities/misc.py b/OpenPNM/Utilities/misc.py
index 5a3889626..0f1724e0c 100644
--- a/OpenPNM/Utilities/misc.py
+++ b/OpenPNM/Utilities/misc.py
@@ -229,9 +229,13 @@ def conduit_lengths(network, throats=None, mode='pore'):
         # Find the pore-to-pore distance, minus the throat length
         lengths = _sp.sqrt(_sp.sum(_sp.square(pcoords[Ps[:, 0]] -
                                    pcoords[Ps[:, 1]]), 1)) - network['throat.length']
-        #   Calculate the fraction of that distance from the first pore
+        lengths[lengths <= 0.0] = 2e-9
+        # Calculate the fraction of that distance from the first pore
         try:
             fractions = pdia[Ps[:, 0]]/(pdia[Ps[:, 0]] + pdia[Ps[:, 1]])
+            # Don't allow zero lengths
+            fractions[fractions == 0.0] = 0.5
+            fractions[fractions == 1.0] = 0.5
         except:
             fractions = 0.5
         plen1 = lengths*fractions
diff --git a/OpenPNM/Utilities/vertexops.py b/OpenPNM/Utilities/vertexops.py
index b7f4a1469..b7b5ce8e5 100644
--- a/OpenPNM/Utilities/vertexops.py
+++ b/OpenPNM/Utilities/vertexops.py
@@ -1,313 +1,18 @@
 import numpy as np
 from scipy.spatial import ConvexHull
 from OpenPNM.Utilities import transformations as tr
-from OpenPNM.Utilities import misc as misc
-from math import atan2
-
-
-def get_throat_geom(verts, normal, fibre_rad):
-    r"""
-    For one set of vertices defining a throat return the key properties
-    This is the main loop for calling other sub-routines.
-    General Method:
-        For each connection or throat defined by the shared vertices
-        Rotate the vertices to align with the xy-plane and get rid of z-coordinate
-        Compute the convex hull of the 2D points giving a set of simplices which
-        define neighbouring vertices in a clockwise fashion.
-        For each triplet calculate the offset position given the fibre radius
-        Check for overlapping vertices and ones that lie outside the original hull -
-        recalculate position to offset from or ignore if all overlapping.
-        Calculate Area and Perimeter if successfully generated offset vertices to
-        replicate eroded throat.
-        Translate back into 3D
-    Any Errors encountered result in the throat area being zero and no vertices being
-    passed back.
-    These Errors are not coding mistakes but failures to obtain an eroded facet with
-    non-zero area:
-    Error 1: Less than 3 vertices in the convex hull - Should never happen
-             (unless 2 points are incredibly close together)
-    Error 2: The largest span of points is less than twice the fibre radius
-             (i.e. throat will definitley be occluded)
-    Error 3: All the offset vertices overlap with at least one other vertex -
-             Throat fully occluded
-    Error 4: Not enough offset vertices to continue - Throat fully occluded
-    Error 5: An offset vertex is outside the original set of points - Throat fully
-             occluded
-    """
-    z_axis = [0, 0, 1]
-    throat_area = 0.0
-    throat_perimeter = 0.0
-    throat_COM = np.zeros([1, 3])
-    output_offset = []
-    Error = 0
-    # For boundaries some facets will already be aligned with the axis -
-    # if this is the case a rotation is unnecessary and could also cause problems
-    angle = tr.angle_between_vectors(normal, z_axis)
-    if angle == 0.0 or angle == np.pi:
-        # We are already aligned
-        rotate_input = False
-        facet = verts
-    else:
-        rotate_input = True
-        M = tr.rotation_matrix(tr.angle_between_vectors(normal, z_axis),
-                               tr.vector_product(normal, z_axis))
-        facet = np.dot(verts, M[:3, :3].T)
-    x = facet[:, 0]
-    y = facet[:, 1]
-    z = facet[:, 2]
-    # Work out span of points and set axes scales to cover this and be equal
-    # in both dimensions
-    x_range = x.max() - x.min()
-    y_range = y.max() - y.min()
-    if (x_range > y_range):
-        my_range = x_range
-    else:
-        my_range = y_range
-    if np.around(z.std(), 3) != 0.000:
-        print('Rotation failed')
-    facet_coords_2D = np.column_stack((x, y))
-    hull = ConvexHull(facet_coords_2D, qhull_options='QJ Pp')
-    verts_2D = facet_coords_2D[hull.vertices]
-    offset = outer_offset(verts_2D, fibre_rad)
-    # At this point we may have overlapping areas for which we need to offset
-    # from a new point
-    overlap_array, sweep_radius, line_points = set_overlap(verts_2D, offset)
-    temp_vert_list = []
-    if len(verts_2D) < 3:
-        # Error: Fused Too Many Verts
-        Error = 1
-    elif(my_range < fibre_rad*2):
-        # Error: Facet Too small to Erode
-        Error = 2
-    else:
-        if all_overlap(overlap_array) is True:
-            # If all overlaps then throat is fully occluded
-            # Error: Throat fully occluded
-            Error = 3
-        else:
-            # If one or two sets of overlaps exist and at least one vertex is not
-            # overlapped then we need to do a bit more work. Do some linalg to find
-            # a new point to offset from saving un-overlapped verts and newly
-            # created verts in a temporary list
-            count = 0
-            temp_verts = verts_2D
-            while True:
-                temp_vert_list = []
-                for i in range(np.shape(line_points)[0]):
-                    if np.sum(overlap_array[i]) == 0.0:
-                        temp_vert_list.append(temp_verts[i])
-                    else:
-                        my_lines = []
-                        for j in range(np.shape(line_points)[0]):
-
-                            if overlap_array[i][j] == 1 and overlap_array[j][i] == 1:
-                                list_a = line_points[i][j]
-                                list_b = line_points[j][i]
-                                my_lines = symmetric_difference(list_a, list_b)
-
-                        my_lines = np.asarray(my_lines)
-
-                        if len(my_lines) == 2:
-                            try:
-                                quad_points = temp_verts[my_lines]
-                                my_new_point = new_point(quad_points)
-                                temp_vert_list.append(my_new_point)
-                            except IndexError:
-                                print('IndexError: ' + str(my_lines))
-                            except TypeError:
-                                print('TypeError: ' + str(my_lines))
-
-                temp_verts = np.asarray(misc.unique_list(temp_vert_list))
-                offset = outer_offset(temp_verts, fibre_rad)
-                overlap_array, sweep_radius, line_points = \
-                    set_overlap(temp_verts, offset)
-
-                if overlap_array.any() is False:
-                    break
-                elif all_overlap(overlap_array) is True:
-                    Error = 3
-                    break
-                elif len(temp_verts) < 3:
-                    Error = 4
-                    break
-                else:
-                    count += 1
-                    temp_verts = np.asarray(fuse_verts(verts=temp_verts,
-                                                       percentage=0.05*count))
-                    offset = outer_offset(temp_verts, fibre_rad)
-                    overlap_array, sweep_radius, line_points = \
-                        set_overlap(temp_verts, offset)
-                    # Continue looping until one of the above conditions is true
-                    # or counter reaches 10
-                if count >= 10:
-                    break
-
-    if len(offset) >= 3 and Error == 0:
-        # Now also check whether any of the offset points lie outside the
-        # original convex hull
-        original_area = np.around(PolyArea2D(verts_2D), 10)
-        all_points = np.concatenate((verts_2D, offset), axis=0)
-        try:
-            total_hull = ConvexHull(all_points, qhull_options='QJ Pp')
-            total_area = np.around(PolyArea2D(all_points[total_hull.vertices]), 10)
-        except np.spatial.qhull.QhullError:
-            print(all_points)
-            total_area = 999
-            Error = 5
-        offset_hull = ConvexHull(offset, qhull_options='QJ Pp')
-        offset_verts_2D = offset[offset_hull.vertices]
-        if total_area > original_area:
-            if Error != 5:
-                Error = 6
-        else:
-            throat_area = PolyArea2D(offset_verts_2D)
-            throat_perimeter = PolyPerimeter2D(offset_verts_2D)
-            throat_COM_2D = PolyWeightedCentroid2D(offset_verts_2D)
-            throat_COM = np.hstack((throat_COM_2D, z[0]))
-        # Make 3D again in rotated plane
-        offset_verts_3D = np.column_stack((offset_verts_2D,
-                                           z[0:len(offset_verts_2D)]))
-        # Get matrix to un-rotate the co-ordinates back to the original orientation
-        # if we rotated in the first place
-        if (rotate_input):
-            M1 = tr.inverse_matrix(M)
-            # Unrotate the offset coordinates
-            output_offset = np.dot(offset_verts_3D, M1[:3, :3].T)
-            throat_COM = np.dot(throat_COM, M1[:3, :3].T)
-        else:
-            output_offset = offset_verts_3D
-
-    return throat_area, throat_perimeter, output_offset, throat_COM, Error
-
-
-def outer_offset(verts, fibre_rad):
-    r"""
-    Routine to loop through all verts and calculate offset position based on
-    neighbours either side. Verts must be in hull order
-    """
-    offset = []
-    for i, vert in enumerate(verts):
-        # Collect three adjacent points and compute the offset of the first
-        triplet = (vert, np.roll(verts, -1, axis=0)[i], np.roll(verts, 1, axis=0)[i])
-        offset.append(offset_vertex(triplet, fibre_rad))
-    offset = np.asarray(offset)
-
-    return offset
-
-
-def offset_vertex(points, rad=0.01):
-    # We are passed in a set of 3 points forming vertices of two adjoining simplexes
-    # of the convex hull of a voronoi facet. We need to offset the vertices normal to
-    # the fibre direction (or adjoining vectors) by the fibre radius. This is
-    # achieved by finding the half angle between the two adjoining vectors and a
-    # direction. Mid-point must be the first in the array
-    p0 = np.array(points[0])
-    p1 = np.array(points[1])
-    p2 = np.array(points[2])
-    # Now make the midpoint the origin
-    vector1 = p1-p0
-    vector2 = p2-p0
-
-    # Find what quadrant the vector is pointing in - atan2 function takes account of
-    # signs. 0 means aligned with x-axis, pi is aligned with -xaxis, positive numbers
-    # are positive y and negative numbers are negative y. The angle between the
-    # vectors should always be within 180 degrees of one another in a convex hull.
-
-    q1 = atan2(vector1[1], vector1[0])
-    q2 = atan2(vector2[1], vector2[0])
-    alpha = 0.5*tr.angle_between_vectors(vector1, vector2)
-
-    # We always want to offset from the first vertex we get to - going anti-clockwise
-    # from the x-axis Check if both vectors point up or both point down - if so the
-    # first one we get to will have smaller q value
-    if q1*q2 >= 0.0:
-        if q1 < q2:
-            theta = q1
-        else:
-            theta = q2
-    else:
-        # If vector 1 is more rotated away from positive xaxis than vector 2 is
-        # rotated away from negative xaxis - use it, and vice-versa
-        if (abs(q1) + abs(q2)) > np.pi:
-            # Vectors are pointing negative x so take whichever has positive q-value
-            # like a pacman facing left
-            if q1 >= 0:
-                theta = q1
-            else:
-                theta = q2
-        else:
-            # Vectors are pointing positive x so take whichever is negative
-            if q1 <= 0:
-                theta = q1
-            else:
-                theta = q2
-
-    # This may cause problems in terms of which way to offset.
-    if alpha == 0:
-        x = 0
-        y = 0
-    else:
-        x = rad*np.cos(alpha+theta)/np.sin(alpha)
-        y = rad*np.sin(alpha+theta)/np.sin(alpha)
-
-    # Add the midpoint back in
-    output = [x+p0[0], y+p0[1]]
-
-    return output
-
-
-def dist2(p1, p2):
-    r"""
-    Pythagoras to compute the square of the distance between two points (in 2D)
-    """
-    return (p1[0] - p2[0])**2 + (p1[1] - p2[1])**2
-
-
-def fuse(points, d):
-    r"""
-    Fuse points together wihin a certain range
-    """
-    ret = []
-    d2 = d * d
-    n = len(points)
-    taken = [False] * n
-    for i in range(n):
-        if not taken[i]:
-            count = 1
-            point = [points[i][0], points[i][1]]
-            taken[i] = True
-            for j in range(i+1, n):
-                if dist2(points[i], points[j]) < d2:
-                    point[0] += points[j][0]
-                    point[1] += points[j][1]
-                    count += 1
-                    taken[j] = True
-            point[0] /= count
-            point[1] /= count
-            ret.append((point[0], point[1]))
-    return ret
-
-
-def fuse_verts(verts, percentage=0.05):
-    r"""
-    Work out the span of the points and therefore the range for fusing them together
-    then call fuse
-    """
-    # Work out largest span
-    x_span = max(verts[:, 0]) - min(verts[:, 0])
-    y_span = max(verts[:, 1]) - min(verts[:, 1])
-    if x_span > y_span:
-        tolerance = x_span*percentage
-    else:
-        tolerance = y_span*percentage
-    # Fuse vertices lying within 5% of the largest span
-    return fuse(verts, tolerance)
 
 
 def PolyArea2D(pts):
     r"""
     returns the area of a 2D polygon given the set of points defining the convex hull
     in correct order
+    Example
+    ---------
+    >>> import OpenPNM.Utilities.vertexops as vo
+    >>> tri = np.array([[0,0],[1,2],[2,0]])
+    >>> vo.PolyArea2D(tri) == 2.0
+    True
     """
     lines = np.hstack([pts, np.roll(pts, -1, axis=0)])
     area = 0.5*abs(sum(x1*y2-x2*y1 for x1, y1, x2, y2 in lines))
@@ -318,6 +23,12 @@ def PolyPerimeter2D(pts):
     r"""
     returns the perimeter of a 2D polygon given the set of points defining the convex
     hull in correct order
+    Example
+    ---------
+    >>> import OpenPNM.Utilities.vertexops as vo
+    >>> quad = np.array([[0,0],[0,1],[1,1],[1,0]])
+    >>> vo.PolyPerimeter2D(quad) == 4.0
+    True
     """
     lines = np.hstack([pts, np.roll(pts, -1, axis=0)])
     perimeter = sum(np.sqrt((x2-x1)**2+(y2-y1)**2) for x1, y1, x2, y2 in lines)
@@ -326,8 +37,14 @@ def PolyPerimeter2D(pts):
 
 def PolyWeightedCentroid2D(pts):
     r"""
-    returns the perimeter of a 2D polygon given the set of points defining the convex
+    returns the centroid of a 2D polygon given the set of points defining the convex
     hull in correct order
+    Example
+    ---------
+    >>> import OpenPNM.Utilities.vertexops as vo
+    >>> quad = np.array([[0,0],[0,2],[2,2],[2,0]])
+    >>> vo.PolyWeightedCentroid2D(quad) == [1.0,1.0]
+    True
     """
     lines = np.hstack([pts, np.roll(pts, -1, axis=0)])
     twice_area = 0.0
@@ -345,193 +62,8 @@ def PolyWeightedCentroid2D(pts):
     return [Cx, Cy]
 
 
-def symmetric_difference(list_a, list_b):
-    r"""
-    Return the combination of two lists without common elements (necessary as sets
-    cannot contain mutable objects)
-    """
-    sym_diff = []
-    sorted_list_a = np.sort(list_a)
-    sorted_list_b = np.sort(list_b)
-    # Add elements in list a if not in list b
-    for element_a in sorted_list_a:
-        match = False
-        for element_b in sorted_list_b:
-            if all(element_a == element_b):
-                match = True
-        if match is False:
-            sym_diff.append(element_a)
-    # Add elements in list b if not in list a
-    for element_b in sorted_list_b:
-        match = False
-        for element_a in sorted_list_a:
-            if all(element_a == element_b):
-                match = True
-        if match is False:
-            sym_diff.append(element_b)
-    return sym_diff
-
-
-def new_point(pairs):
-    r"""
-    Passed 2 pairs of points defining lines either side of overlapped offset vertices
-    need to calculate the new point to offset from given the orientation of the outer
-    fibres
-    """
-    m1, c1 = line_equation(pairs[0])
-    m2, c2 = line_equation(pairs[1])
-
-    if m1 == np.inf:
-        # line1 is a straight line x = c1
-        x = c1
-        y = (m2*c1) + c2
-    elif (m2 == np.inf):
-        # line2 is a straight line x = c2
-        x = c2
-        y = (m1*c2) + c1
-    else:
-        try:
-            x = (c2-c1) / (m1-m2)
-            y = (m1*c2 - m2*c1) / (m1-m2)
-        except RuntimeWarning:
-            x = 0
-            y = 0
-    return x, y
-
-
-def line_equation(points):
-    r"""
-    Return the gradient and y intercept of a straight line given 2 points
-    """
-    x_coords, y_coords = zip(*points)
-    dy = y_coords[1]-y_coords[0]
-    dx = x_coords[1]-x_coords[0]
-    if dx == 0:
-        m = np.inf
-        c = x_coords[1]
-    else:
-        m = dy / dx
-        c = y_coords[1] - m*x_coords[1]
-
-    return m, c
-
-
-def set_overlap(verts, offset):
-    r"""
-    Given a set of vertices and a set of offset vertices, evaluate whether any of the
-    offset vertices overlap. This is then used to recalculate points from which to
-    offset
-    """
-    dim = len(verts)
-    overlap_array = np.zeros(shape=(dim, dim))
-    sweep_radius = np.zeros(len(verts))
-    for i, zone_centre in enumerate(verts):
-        sweep_radius[i] = np.sqrt(dist2(zone_centre, offset[i]))
-        for j, test_point in enumerate(offset):
-            test_radius = np.sqrt(dist2(zone_centre, test_point))
-            if (test_radius < sweep_radius[i]):
-                # Fill in both so that they are both recalculated later i overlapping
-                # j doesn't necessarily mean j overlaps i
-                overlap_array[i][j] = 1
-                overlap_array[j][i] = 1
-    # Join up overlapping regions of points
-    for i in range(dim):
-        for j in range(dim):
-            # If an overlap exist look at what others exist for that vertex
-            if overlap_array[i][j] == 1:
-                for k in range(dim):
-                    if overlap_array[j][k] == 1 and k != i:
-                        overlap_array[i][k] = 1
-
-    line_points_out = line_points(overlap_array)
-
-    return overlap_array, sweep_radius, line_points_out
-
-
-def line_points(array):
-    r"""
-    We are passed a square array containing a list of overlap results. rows represent
-    vertices and columns represent offset vertices. If an overlap occurs in the span
-    between offset j and offset i then a 1 will result in [i][j]. As the vertices are
-    in hull order and our aim is to create a new point from which to offset using
-    connected fibres we want to identify the correct points to use to define our
-    lines.
-
-      e__________ d
-       |        | c
-       |       /
-       |     /
-       |   /
-       |_/
-       a b
-    If we have 5 points in a hull a,b,c,d,e where a overlaps with b and c overlaps
-    with d (and visa versa) the array will look like:
-    [0,1,0,0,0]
-    [1,0,0,0,0]
-    [0,0,0,1,0]
-    [0,0,1,0,0]
-    [0,0,0,0,0]
-
-    This means that c and d are within a fibre's width of each other and the line
-    between them does not represent the fibre. Instead we want to extend the outer
-    lines (bc and de) to see where they would meet and offset from this point. Roll
-    up and down to find the first unoverlapped index from which to start each line
-    from then go back one to get the two points to form a line.
-    """
-    dim = np.shape(array)[0]
-    index = range(dim)
-    line_points = np.ndarray(shape=[dim, dim], dtype=object)
-    for i in range(dim):
-        for j in range(dim):
-            if array[i][j] == 1 and array[j][i] == 1:
-                # Roll forwards to find the first unoverlapped index
-                k = 1
-                while k < dim:
-                    if np.roll(array[i], -k, axis=0)[j] == 0:
-                        break
-                    else:
-                        k += 1
-                # Save the indices of the first unoverlapped index and the one before
-                # to create the line
-                forward_line = [np.roll(index, -k, axis=0)[j],
-                                np.roll(index, -(k-1), axis=0)[j]]
-                forward_line.sort()
-                # Roll backwards to find the first unoverlapped index
-                k = 1
-                while k < dim:
-                    if np.roll(array[i], k, axis=0)[j] == 0:
-                        break
-                    else:
-                        k += 1
-                # Save the indices of the first unoverlapped index and the one before
-                # to create the line
-                backward_line = [np.roll(index, k, axis=0)[j],
-                                 np.roll(index, (k - 1), axis=0)[j]]
-                backward_line.sort()
-                line_points[i][j] = (forward_line, backward_line)
-
-    return line_points
-
-
-def all_overlap(array):
-    r"""
-    Find out whether all offset vertices (columns) are overlapped by at least one
-    other. If so then throat is fully occluded
-    """
-    dim = np.shape(array)[0]
-    overlap = [False]*dim
-    all_overlap = False
-    for i in range(dim):
-        for j in range(dim):
-            if array[j][i] == 1:
-                overlap[i] = True
-    if sum(overlap) == dim:
-        all_overlap = True
-
-    return all_overlap
-
-
-def scale(network, scale_factor=[1, 1, 1], preserve_vol=True):
+def scale(network, scale_factor=[1, 1, 1], preserve_vol=False,
+          linear_scaling=[False, False, False]):
     r"""
     A method for scaling the coordinates and vertices to create anisotropic networks
     The original domain volume can be preserved by setting preserve_vol = True
@@ -546,32 +78,79 @@ def scale(network, scale_factor=[1, 1, 1], preserve_vol=True):
     >>> B1 = pn.pores("left_boundary")
     >>> B2 = pn.pores("right_boundary")
     >>> Vol = vo.vertex_dimension(pn,B1,B2)
-    >>> vo.scale(network=pn,scale_factor=[2,1,1])
+    >>> vo.scale(network=pn,scale_factor=[2,1,1],preserve_vol=True)
     >>> Vol2 = vo.vertex_dimension(pn,B1,B2)
     >>> np.around(Vol-Vol2,5) == 0.0
     True
     >>> vo.scale(network=pn,scale_factor=[2,1,1],preserve_vol=False)
     >>> Vol3 = vo.vertex_dimension(pn,B1,B2)
-    >>> np.around(Vol3/Vol,5)
-    2.0
+    >>> np.around(Vol3/Vol,5) == 2.0
+    True
 
     """
     from scipy.special import cbrt
     import scipy as sp
+    minmax = np.around(vertex_dimension(network=network,
+                                        face1=network.pores(), parm='minmax'), 10)
     scale_factor = np.asarray(scale_factor)
     if preserve_vol is True:
         scale_factor = scale_factor/(cbrt(sp.prod(scale_factor)))
-    network['pore.coords'] = network['pore.coords']*scale_factor
+
+    lin_scale = _linear_scale_factor(network["pore.coords"], minmax,
+                                     scale_factor, linear_scaling)
+
+    network["pore.coords"] = network["pore.coords"]*lin_scale
     # Cycle through all vertices of all pores updating vertex values
     for pore in network.pores():
         for i, vert in network['pore.vert_index'][pore].items():
-            network['pore.vert_index'][pore][i] = \
-                network['pore.vert_index'][pore][i]*scale_factor
+            vert_scale = _linear_scale_factor(vert, minmax,
+                                              scale_factor, linear_scaling)
+            network["pore.vert_index"][pore][i] = vert*vert_scale
     # Cycle through all vertices of all throats updating vertex values
     for throat in network.throats():
         for i, vert in network['throat.vert_index'][throat].items():
-            network['throat.vert_index'][throat][i] = \
-                network['throat.vert_index'][throat][i]*scale_factor
+
+            vert_scale = _linear_scale_factor(vert, minmax,
+                                              scale_factor, linear_scaling)
+            network["throat.vert_index"][throat][i] = vert*vert_scale
+    # Scale the vertices on the voronoi diagram stored on the network
+    # These are used for adding boundaries on the Delaunay network class
+    vert = network._vor.vertices
+    vert_scale = _linear_scale_factor(vert, minmax, scale_factor, linear_scaling)
+    network._vor.vertices = vert*vert_scale
+
+
+def _linear_scale_factor(points=None,
+                         minmax=[0, 1, 0, 1, 0, 1],
+                         scale_factor=[1, 1, 1],
+                         linear_scaling=[False, False, False]):
+    r"""
+    Work out the linear scale factor of a point or set of points based on the
+    domain extent, an absolute scale factor and linear_scaling booleans for each
+    coordinate. If all False the absolute scaling is applied equally across the
+    domain. If one linear_scaling boolean is True then a linear function is
+    applied to scaling the co-ordinates along that axis. If more than one boolean
+    is true then a combined linear function is applied
+    """
+    [xmin, xmax, ymin, ymax, zmin, zmax] = minmax
+    max_array = np.array([(xmax-xmin), (ymax-ymin), (zmax-zmin)])
+    pos_array = points/max_array
+    shape = np.shape(points)
+    if len(shape) == 1:
+        combined_pos = np.ones([1])
+        for i in range(3):
+            if linear_scaling[i] == True:
+                combined_pos *= pos_array[i]
+        lin_scale = (scale_factor - 1)*combined_pos + 1
+    else:
+        combined_pos = np.ones([shape[0]])
+        for i in range(3):
+            if linear_scaling[i] == True:
+                combined_pos *= pos_array[:, i]
+        pos_array = np.vstack((combined_pos, combined_pos, combined_pos)).T
+        lin_scale = (scale_factor - 1)*pos_array + 1
+
+    return lin_scale
 
 
 def vertex_dimension(network, face1=[], face2=[], parm='volume'):
@@ -682,11 +261,11 @@ def porosity(network):
     """
     domain_vol = vertex_dimension(network, network.pores(), parm='volume')
     try:
-        pore_vol = sum(network['pore.volume'])
+        pore_vol = np.sum(network['pore.volume'])
     except KeyError:
         print('Geometries must be assigned first')
         pore_vol = 0
-    porosity = pore_vol/domain_vol
+    porosity = np.around(pore_vol/domain_vol, 3)
     return porosity
 
 
@@ -735,29 +314,32 @@ def tortuosity(network=None):
     return 1 / np.cos(np.array([theta_x, theta_y, theta_z]))
 
 
-def print_throat(geom, throats_in):
+def print_throat(geometry, throats, fig=None):
     r"""
     Print a given throat or list of throats accepted as [1, 2, 3, ..., n]
-    e.g geom.print_throat([34, 65, 99])
     Original vertices plus offset vertices are rotated to align with
     the z-axis and then printed in 2D
+    e.g vo.print_throat(geom, [34, 65, 99])
     """
     import matplotlib.pyplot as plt
-    throats = []
-    for throat in throats_in:
-        if throat in range(geom.num_throats()):
-            throats.append(throat)
+    throat_list = []
+    for throat in throats:
+        if throat in range(geometry.num_throats()):
+            throat_list.append(throat)
         else:
             print('Throat: ' + str(throat) + ' not part of geometry')
-    if len(throats) > 0:
-        verts = geom['throat.vertices'][throats]
-        offsets = geom['throat.offset_vertices'][throats]
-        normals = geom['throat.normal'][throats]
-        coms = geom['throat.centroid'][throats]
-        incentre = geom['throat.incentre'][throats]
-        inradius = 0.5*geom['throat.indiameter'][throats]
-        for i in range(len(throats)):
-            fig = plt.figure()
+    if len(throat_list) > 0:
+        verts = geometry['throat.vertices'][throat_list]
+        offsets = geometry['throat.offset_vertices'][throat_list]
+        normals = geometry['throat.normal'][throat_list]
+        coms = geometry['throat.centroid'][throat_list]
+        incentre = geometry['throat.incentre'][throat_list]
+        inradius = 0.5*geometry['throat.indiameter'][throat_list]
+        row_col = np.ceil(np.sqrt(len(throat_list)))
+        for i in range(len(throat_list)):
+            if fig is None:
+                fig = plt.figure()
+            ax = fig.add_subplot(row_col, row_col, i+1)
             vert_2D = tr.rotate_and_chop(verts[i], normals[i], [0, 0, 1])
             hull = ConvexHull(vert_2D, qhull_options='QJ Pp')
             for simplex in hull.simplices:
@@ -795,43 +377,47 @@ def print_throat(geom, throats_in):
             u = inradius[i]*np.cos(t)+incent[0][0]
             v = inradius[i]*np.sin(t)+incent[0][1]
             plt.plot(u, v, 'r-')
+            ax.ticklabel_format(style='sci', scilimits=(0, 0))
             fig.show()
     else:
-        print('Please provide throat indices')
+        print("Please provide throat indices")
+    return fig
 
 
-def print_pore(geom, pores, fig=None, axis_bounds=None):
+def print_pore(geometry, pores, fig=None, axis_bounds=None, include_points=False):
     r"""
     Print all throats around a given pore or list of pores accepted
     as [1, 2, 3, ..., n]
-    e.g geom.print_pore([34, 65, 99])
+    e.g vo.print_pore(geom, [34, 65, 99])
     Original vertices plus offset vertices used to create faces and
     then printed in 3D
     To print all pores (n)
     pore_range = np.arange(0,n-1,1)
-    geom.print_pore(pore_range)
+    vo.print_pore(geom, pore_range)
     """
     import matplotlib.pyplot as plt
     from mpl_toolkits.mplot3d import Axes3D
     from mpl_toolkits.mplot3d.art3d import Poly3DCollection
-    return_fig = False
     if len(pores) > 0:
-        net_pores = geom.map_pores(geom._net, pores)
-        centroids = geom['pore.centroid'][pores]
-        throats = geom._net.map_throats(geom,
-                                        net_throats,
-                                        return_mapping=True)['target']
+        net_pores = geometry.map_pores(geometry._net, pores)
+        centroids = geometry['pore.centroid'][pores]
+        coords = geometry._net['pore.coords'][net_pores]
+        net_throats = geometry._net.find_neighbor_throats(pores=net_pores)
+        throats = geometry._net.map_throats(geometry,
+                                            net_throats,
+                                            return_mapping=True)['target']
+        tcentroids = geometry["throat.centroid"][throats]
         # Can't create volume from one throat
         if 1 <= len(throats):
-            verts = geom['throat.vertices'][throats]
-            normals = geom['throat.normal'][throats]
+            verts = geometry['throat.vertices'][throats]
+            normals = geometry['throat.normal'][throats]
             # Get verts in hull order
             ordered_verts = []
             for i in range(len(verts)):
                 vert_2D = tr.rotate_and_chop(verts[i], normals[i], [0, 0, 1])
                 hull = ConvexHull(vert_2D, qhull_options='QJ Pp')
                 ordered_verts.append(verts[i][hull.vertices])
-            offsets = geom['throat.offset_vertices'][throats]
+            offsets = geometry['throat.offset_vertices'][throats]
             ordered_offs = []
             for i in range(len(offsets)):
                 offs_2D = tr.rotate_and_chop(offsets[i], normals[i], [0, 0, 1])
@@ -840,13 +426,11 @@ def print_pore(geom, pores, fig=None, axis_bounds=None):
             # Get domain extents for setting axis
             if axis_bounds is None:
                 [xmin, xmax, ymin, ymax, zmin, zmax] = \
-                    vertex_dimension(geom._net, pores, parm='minmax')
+                    vertex_dimension(geometry._net, pores, parm='minmax')
             else:
                 [xmin, xmax, ymin, ymax, zmin, zmax] = axis_bounds
             if fig is None:
                 fig = plt.figure()
-            else:
-                return_fig is True
             ax = fig.gca(projection='3d')
             outer_items = Poly3DCollection(ordered_verts, linewidths=1,
                                            alpha=0.2, zsort='min')
@@ -861,11 +445,15 @@ def print_pore(geom, pores, fig=None, axis_bounds=None):
             ax.set_xlim(xmin, xmax)
             ax.set_ylim(ymin, ymax)
             ax.set_zlim(zmin, zmax)
-            ax.scatter(centroids[:, 0], centroids[:, 1], centroids[:, 2], c='y')
+            if include_points:
+                ax.scatter(centroids[:, 0], centroids[:, 1], centroids[:, 2], c='y')
+                ax.scatter(tcentroids[:, 0], tcentroids[:, 1], tcentroids[:, 2],
+                           c='r')
+                ax.scatter(coords[:, 0], coords[:, 1], coords[:, 2], c='b')
+            ax.ticklabel_format(style='sci', scilimits=(0, 0))
             plt.show()
         else:
             print_throat(throats)
     else:
         print('Please provide pore indices')
-    if return_fig is True:
-        return fig
+    return fig
diff --git a/article recreation/gostick_recreation.py b/article recreation/gostick_recreation.py
index a9d4ba724..4c30aa851 100644
--- a/article recreation/gostick_recreation.py
+++ b/article recreation/gostick_recreation.py
@@ -4,7 +4,7 @@ import matplotlib.pyplot as plt
 Lc = 40.5e-6
 
 #1 setting up network
-sgl = OpenPNM.Network.Cubic([26, 26, 10], spacing=Lc, name='SGL10BA')
+sgl = OpenPNM.Network.Cubic([10, 10, 10], spacing=Lc, name='SGL10BA')
 sgl.add_boundaries()
 
 #2 set up geometries
@@ -111,17 +111,17 @@ for x in range(num_seq+1):
 
         #run Stokes Flow and find Permeability
         #single phase
-        Stokes_alg_single_phase_air = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_single_phase_air',network=sgl,phase=air)
-        Stokes_alg_single_phase_water = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_single_phase_water',network=sgl,phase=water)
+        Stokes_alg_single_phase_air = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_single_phase_air'+str(bound_increment)+str(x),network=sgl,phase=air)
+        Stokes_alg_single_phase_water = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_single_phase_water'+str(bound_increment)+str(x),network=sgl,phase=water)
         
-        Fickian_alg_single_phase_air = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_single_phase_air',network=sgl,phase=air)
-        Fickian_alg_single_phase_water = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_single_phase_water',network=sgl,phase=water)
+        Fickian_alg_single_phase_air = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_single_phase_air'+str(bound_increment)+str(x),network=sgl,phase=air)
+        Fickian_alg_single_phase_water = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_single_phase_water'+str(bound_increment)+str(x),network=sgl,phase=water)
         
-        Stokes_alg_multi_phase_air = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_multi_phase_air',network=sgl,phase=air)
-        Stokes_alg_multi_phase_water = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_multi_phase_water',network=sgl,phase=water)
+        Stokes_alg_multi_phase_air = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_multi_phase_air'+str(bound_increment)+str(x),network=sgl,phase=air)
+        Stokes_alg_multi_phase_water = OpenPNM.Algorithms.StokesFlow(name='Stokes_alg_multi_phase_water'+str(bound_increment)+str(x),network=sgl,phase=water)
         
-        Fickian_alg_multi_phase_air = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_multi_phase_air',network=sgl,phase=air)
-        Fickian_alg_multi_phase_water = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_multi_phase_water',network=sgl,phase=water)
+        Fickian_alg_multi_phase_air = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_multi_phase_air'+str(bound_increment)+str(x),network=sgl,phase=air)
+        Fickian_alg_multi_phase_water = OpenPNM.Algorithms.FickianDiffusion(name='Fickian_alg_multi_phase_water'+str(bound_increment)+str(x),network=sgl,phase=water)
         
         BC1_pores = sgl.pores(labels=bounds[bound_increment][0]+'_boundary')
         BC2_pores = sgl.pores(labels=bounds[bound_increment][1]+'_boundary')
diff --git a/run_script.py b/run_script.py
index 93f931535..39fd15319 100644
--- a/run_script.py
+++ b/run_script.py
@@ -1,6 +1,8 @@
 import OpenPNM
+import time
+st = time.time()
 print('-----> Using OpenPNM version: '+OpenPNM.__version__)
-
+from OpenPNM.Geometry import models as gm
 ctrl = OpenPNM.Base.Controller()
 #==============================================================================
 '''Build Topological Network'''
@@ -14,7 +16,7 @@ pn.add_boundaries()
 Ps = pn.pores('boundary',mode='not')
 Ts = pn.find_neighbor_throats(pores=Ps,mode='intersection',flatten=True)
 geom = OpenPNM.Geometry.Toray090(network=pn,pores=Ps,throats=Ts)
-
+geom.models.add(propname='throat.length',model=gm.throat_length.straight)
 Ps = pn.pores('boundary')
 Ts = pn.find_neighbor_throats(pores=Ps,mode='not_intersection')
 boun = OpenPNM.Geometry.Boundary(network=pn,pores=Ps,throats=Ts)
@@ -101,3 +103,4 @@ except Exception as e:
 #------------------------------------------------------------------------------
 ctrl.export()
 
+print("sim time:" + str(time.time()-st))
diff --git a/setup.py b/setup.py
index 5daa826c5..fafb005bc 100644
--- a/setup.py
+++ b/setup.py
@@ -40,7 +40,8 @@ setup(
     install_requires = [
         'numpy',
         'scipy>=0.14.0',
-        'matplotlib'
+        'matplotlib',
+        'scikit-image'
     ],
     author = 'OpenPNM Team',
     author_email = 'jeff.gostick@mcgill.ca',


Test patch:
diff --git a/test/unit/Geometry/VoronoiTest.py b/test/unit/Geometry/VoronoiTest.py
index 3d56c73aa..f90dffbb1 100644
--- a/test/unit/Geometry/VoronoiTest.py
+++ b/test/unit/Geometry/VoronoiTest.py
@@ -1,3 +1,84 @@
+import OpenPNM
+import scipy as sp
+import matplotlib.pyplot as plt
+from OpenPNM.Geometry import models as gm
+import os
+
+
 class VoronoiTest:
-    def test_generate(self):
-        pass
+
+    def setup_class(self):
+        bp = sp.array([[0.2, 0.2, 0.2], [0.2, 0.8, 0.2], [0.8, 0.2, 0.2],
+                       [0.8, 0.8, 0.2], [0.2, 0.2, 0.8], [0.2, 0.8, 0.8],
+                       [0.8, 0.2, 0.8], [0.8, 0.8, 0.8]])
+        scale = 1e-4
+        p = (sp.random.random([len(bp), 3])-0.5)/10000
+        bp += p
+        self.ctrl = OpenPNM.Base.Controller()
+        self.net = OpenPNM.Network.Delaunay(domain_size=[scale, scale, scale],
+                                            base_points=bp*scale)
+        self.net.add_boundaries()
+        Ps = self.net.pores()
+        Ts = self.net.throats()
+        fibre_rad = 5e-6
+        self.geo_vox = OpenPNM.Geometry.Voronoi(network=self.net,
+                                                pores=Ps,
+                                                throats=Ts,
+                                                fibre_rad=fibre_rad,
+                                                voxel_vol=True,
+                                                name='vor_vox')
+
+    def test_props_all(self):
+        a = self.geo_vox.props()
+        assert sorted(a) == ['pore.area', 'pore.centroid', 'pore.diameter',
+                             'pore.fibre_voxels', 'pore.indiameter',
+                             'pore.pore_voxels', 'pore.seed',
+                             'pore.vertices', 'pore.volume',
+                             'throat.area', 'throat.c2c', 'throat.centroid',
+                             'throat.diameter', 'throat.incentre',
+                             'throat.indiameter', 'throat.length',
+                             'throat.normal', 'throat.offset_vertices',
+                             'throat.perimeter', 'throat.seed',
+                             'throat.shape_factor', 'throat.surface_area',
+                             'throat.vertices', 'throat.volume']
+
+    def test_get_fibre_slice(self):
+        slc = self.geo_vox.get_fibre_slice(index=[0, 50, 0])
+        assert sp.shape(slc) == (101, 101)
+
+    def test_compress_geom(self):
+        b1 = self.net.pores('bottom_boundary')
+        b2 = self.net.pores('top_boundary')
+        height1 = self.net.domain_length(b1, b2)
+        self.geo_vox.compress_geometry(factor=[1, 1, 0.5])
+        height2 = self.net.domain_length(b1, b2)
+        assert sp.around(height1/height2, 5) == 2.0
+
+    def test_make_fibre_image(self):
+        del(self.geo_vox._fibre_image)
+        self.geo_vox.make_fibre_image()
+        assert hasattr(self.geo_vox, '_fibre_image') == True
+
+    def test_voronoi_vert(self):
+        self.ctrl.purge_object(self.geo_vox)
+        fibre_rad = 5e-6
+        Ps = self.net.pores()
+        Ts = self.net.throats()
+        self.geo_vert = OpenPNM.Geometry.Voronoi(network=self.net,
+                                                 pores=Ps,
+                                                 throats=Ts,
+                                                 fibre_rad=fibre_rad,
+                                                 voxel_vol=False,
+                                                 name='vor_vert')
+        self.geo_vert.models.add(propname='throat.area2',
+                                 model=gm.throat_area.voronoi)
+        assert sp.all(sp.absolute(1 - self.geo_vert["throat.area2"] /
+                                  self.geo_vert['throat.area']) < 0.05)
+        self.geo_vert.models.add(propname='throat.centroid2',
+                                 model=gm.throat_centroid.centre_of_mass)
+        assert sp.absolute(sp.sum(self.geo_vert['throat.centroid'] -
+                                  self.geo_vert['throat.centroid2'])) < 1e-10
+        self.geo_vert.models.add(propname='throat.perimeter2',
+                                 model=gm.throat_perimeter.voronoi)
+        assert sp.all(sp.absolute(1 - self.geo_vert['throat.perimeter2'] /
+                                  self.geo_vert['throat.perimeter']) < 0.05)
diff --git a/test/unit/Geometry/models/ThroatDiameterTest.py b/test/unit/Geometry/models/ThroatDiameterTest.py
index c5e4e5d1b..e8264cea6 100644
--- a/test/unit/Geometry/models/ThroatDiameterTest.py
+++ b/test/unit/Geometry/models/ThroatDiameterTest.py
@@ -1,3 +1,29 @@
+import OpenPNM
+from OpenPNM.Geometry import models as gm
+import numpy as np
+
+
 class ThroatDiameterTest:
-    def test_cylinder(self):
-        pass
+
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[2, 2, 2])
+        self.net.add_boundaries()
+        Ps = self.net.pores('boundary', mode='not')
+        Ts = self.net.throats('boundary', mode='not')
+        self.geo = OpenPNM.Geometry.TestGeometry(network=self.net,
+                                                 pores=Ps,
+                                                 throats=Ts)
+        BPs = self.net.pores('boundary')
+        BTs = self.net.throats('boundary')
+        self.boun = OpenPNM.Geometry.Boundary(network=self.net,
+                                              pores=BPs,
+                                              throats=BTs)
+
+    def test_min_pore(self):
+        self.geo.models.add(propname='throat.diameter',
+                            model=gm.throat_diameter.minpore)
+        self.boun.models.add(propname='throat.diameter',
+                             model=gm.throat_diameter.minpore)
+        BPs = self.net["throat.conns"][self.net.throats('boundary')][:, 0]
+        assert np.sum(self.net["pore.diameter"][BPs] -
+                      self.boun["throat.diameter"]) == 0.0
diff --git a/test/unit/Network/DelaunayCubicTest.py b/test/unit/Network/DelaunayCubicTest.py
index 4a576b53c..f455d7364 100644
--- a/test/unit/Network/DelaunayCubicTest.py
+++ b/test/unit/Network/DelaunayCubicTest.py
@@ -1,3 +1,41 @@
+import OpenPNM
+import OpenPNM.Utilities.vertexops as vo
+
+
 class DelaunayCubicTest:
-    def test_generate_pores(self):
-        pass
+
+    def test_simple_cubic(self):
+        pn = OpenPNM.Network.DelaunayCubic(shape=[3, 3, 3])
+        pn.add_boundaries()
+        B1 = pn.pores("top_boundary")
+        B2 = pn.pores("bottom_boundary")
+        assert vo.vertex_dimension(pn, B1, B2, 'minmax') == [0.0, 3.0, 0.0,
+                                                             3.0, 0.0, 3.0]
+        assert pn.num_pores() == 81
+
+    def test_orthorhombic(self):
+        pn = OpenPNM.Network.DelaunayCubic(shape=[3, 3, 3], arrangement='O')
+        pn.add_boundaries()
+        B1 = pn.pores("top_boundary")
+        B2 = pn.pores("bottom_boundary")
+        assert vo.vertex_dimension(pn, B1, B2, 'minmax') == [0.0, 3.0, 0.0,
+                                                             3.0, 0.0, 3.0]
+        assert pn.num_pores() == 81
+
+    def test_bcc(self):
+        pn = OpenPNM.Network.DelaunayCubic(shape=[3, 3, 3], arrangement='BCC')
+        pn.add_boundaries()
+        B1 = pn.pores("top_boundary")
+        B2 = pn.pores("bottom_boundary")
+        assert vo.vertex_dimension(pn, B1, B2, 'minmax') == [0.0, 3.0, 0.0,
+                                                             3.0, 0.0, 3.0]
+        assert pn.num_pores() == 89
+
+    def test_fcc(self):
+        pn = OpenPNM.Network.DelaunayCubic(shape=[3, 3, 3], arrangement='FCC')
+        pn.add_boundaries()
+        B1 = pn.pores("top_boundary")
+        B2 = pn.pores("bottom_boundary")
+        assert vo.vertex_dimension(pn, B1, B2, 'minmax') == [0.0, 3.0, 0.0,
+                                                             3.0, 0.0, 3.0]
+        assert pn.num_pores() == 141
diff --git a/test/unit/Network/DelaunayTest.py b/test/unit/Network/DelaunayTest.py
index fac3abb4e..371a693a7 100644
--- a/test/unit/Network/DelaunayTest.py
+++ b/test/unit/Network/DelaunayTest.py
@@ -1,3 +1,41 @@
+import scipy as sp
+import OpenPNM
+import os
+
+
 class DelaunayTest:
-    def test_generate(self):
-        pass
+
+    def setup_class(self):
+        bp = sp.array([[0.2, 0.2, 0.2], [0.2, 0.8, 0.2], [0.8, 0.2, 0.2],
+                       [0.8, 0.8, 0.2], [0.2, 0.2, 0.8], [0.2, 0.8, 0.8],
+                       [0.8, 0.2, 0.8], [0.8, 0.8, 0.8]])
+        scale = 1e-4
+        self.scale = scale
+        p = (sp.random.random([len(bp), 3])-0.5)/10000
+        bp += p
+        self.ctrl = OpenPNM.Base.Controller()
+        self.net = OpenPNM.Network.Delaunay(domain_size=[scale, scale, scale],
+                                            base_points=bp*scale)
+        self.net.add_boundaries()
+
+    def test_trim_occluded_throats(self):
+        Np = self.net.num_pores()
+        Nt = self.net.num_throats()
+        T = self.net.throats('boundary')[0]
+        self.net["throat.area"] = 1
+        self.net["throat.area"][T] = 0
+        self.net.trim_occluded_throats()
+        assert self.net.num_pores() == Np - 1
+        assert self.net.num_throats() == Nt - 1
+
+    def test_domain_dimensions(self):
+        face1 = self.net.pores('bottom_boundary')
+        face2 = self.net.pores('top_boundary')
+        assert sp.around(self.net.domain_length(face1, face2) / self.scale,
+                         3) == 1.0
+        assert sp.around(self.net.domain_area(face1) / self.scale**2,
+                         3) == 1.0
+
+    def test_export_vor_fibres(self):
+        self.net.export_vor_fibres()
+        os.remove('fibres.p')
diff --git a/test/unit/Utilities/VertexOpsTest.py b/test/unit/Utilities/VertexOpsTest.py
index 73ae451da..4aa8754ec 100644
--- a/test/unit/Utilities/VertexOpsTest.py
+++ b/test/unit/Utilities/VertexOpsTest.py
@@ -1,3 +1,50 @@
+import OpenPNM
+import scipy as sp
+import matplotlib.pyplot as plt
+import OpenPNM.Utilities.vertexops as vo
+
+
 class VertexOpsTest:
-    def test_get_throat_geom(self):
-        pass
+    def setup_class(self):
+        bp = sp.array([[0.2, 0.2, 0.2], [0.2, 0.8, 0.2], [0.8, 0.2, 0.2],
+                       [0.8, 0.8, 0.2], [0.2, 0.2, 0.8], [0.2, 0.8, 0.8],
+                       [0.8, 0.2, 0.8], [0.8, 0.8, 0.8]])
+        scale = 1e-4
+        p = (sp.random.random([len(bp), 3])-0.5)/10000
+        bp += p
+        self.ctrl = OpenPNM.Base.Controller()
+        self.net = OpenPNM.Network.Delaunay(domain_size=[scale, scale, scale],
+                                            base_points=bp*scale)
+        self.net.add_boundaries()
+        Ps = self.net.pores()
+        Ts = self.net.throats()
+        self.fibre_rad = 5e-6
+        self.geo = OpenPNM.Geometry.Voronoi(network=self.net,
+                                            pores=Ps,
+                                            throats=Ts,
+                                            fibre_rad=self.fibre_rad,
+                                            voxel_vol=False,
+                                            name='vor')
+
+    def test_scale(self):
+        factor = [1, 1, 0.5]
+        vo.scale(network=self.net,
+                 scale_factor=factor,
+                 linear_scaling=[True, False, False],
+                 preserve_vol=False)
+
+    def test_porosity(self):
+        por = vo.porosity(self.net)
+        assert por < 1.0
+
+    def test_tortuosity(self):
+        tor = vo.tortuosity(self.net)
+        assert sp.all(tor > 1.0)
+
+    def test_pore2centroid(self):
+        temp_coords = self.net['pore.coords']
+        self.geo['pore.centroid'] = sp.ones([self.geo.num_pores(), 3])
+        vo.pore2centroid(self.net)
+        assert sp.sum(self.net['pore.coords'] -
+                      sp.ones([self.geo.num_pores(), 3])) == 0.0
+        self.net['pore.coords'] = temp_coords