Repository: PMEAL/OpenPNM
Pull request #410 (2015-06-26 15:45:33)

Feature/voxel voronoi
Opening PR for @TomTranter to get his changes on to master


Patch:
diff --git a/OpenPNM/Algorithms/__OrdinaryPercolation__.py b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
index d34e782eb..45c9650f4 100644
--- a/OpenPNM/Algorithms/__OrdinaryPercolation__.py
+++ b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
@@ -258,6 +258,23 @@ class OrdinaryPercolation(GenericAlgorithm):
         self['pore.inv_Pc'][self['pore.trapped'] > 0] = sp.inf
         self['throat.inv_Pc'][self['throat.trapped'] > 0] = sp.inf
 
+    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['pore.inv_Pc']/Pc)**eta
+        Swp[self['pore.inv_Pc'] > 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
@@ -349,6 +366,7 @@ class OrdinaryPercolation(GenericAlgorithm):
         if sp.mean(self._inv_phase['pore.contact_angle']) < 90:
             Snwp_p = 1 - Snwp_p
             Snwp_t = 1 - Snwp_t
+            Snwp_all = 1 - Snwp_all
             PcPoints *= -1
         fig = plt.figure()
         plt.plot(PcPoints, Snwp_all, 'g.-')
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 58e24980d..eb03920d4 100644
--- a/OpenPNM/Geometry/__GenericGeometry__.py
+++ b/OpenPNM/Geometry/__GenericGeometry__.py
@@ -5,10 +5,9 @@ GenericGeometry -- Base class to manage pore scale geometry
 ===============================================================================
 
 """
-
 import scipy as sp
-import OpenPNM.Geometry.models
 from OpenPNM.Base import Core
+from OpenPNM.Postprocessing import Plots
 from OpenPNM.Base import logging
 from OpenPNM.Network import GenericNetwork
 logger = logging.getLogger(__name__)
@@ -33,6 +32,7 @@ class GenericGeometry(Core):
 
     Examples
     --------
+    >>> import OpenPNM
     >>> pn = OpenPNM.Network.TestNet()
     >>> Ps = pn.pores()  # Get all pores
     >>> Ts = pn.throats()  # Get all throats
@@ -69,7 +69,7 @@ class GenericGeometry(Core):
             return super(GenericGeometry, self).__getitem__(key)
         if key == 'throat.conns':  # Handle specifically
             [P1, P2] = \
-                self._net['throat.conns'][self._net[element + '.' + self.name]].T
+                self._net['throat.conns'][self._net[element+'.'+self.name]].T
             Pmap = sp.zeros((self._net.Np,), dtype=int) - 1
             Pmap[self._net.pores(self.name)] = self.Ps
             conns = sp.array([Pmap[P1], Pmap[P2]]).T
@@ -97,9 +97,21 @@ class GenericGeometry(Core):
         """
         pores = self._parse_locations(pores)
         throats = self._parse_locations(throats)
-        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)
+
+    def plot_histograms(self,
+                        throat_diameter='throat.diameter',
+                        pore_diameter='pore.diameter',
+                        throat_length='throat.length'):
+
+        Plots.distributions(obj=self,
+                            throat_diameter=throat_diameter,
+                            pore_diameter=pore_diameter,
+                            throat_length=throat_length)
+
+    plot_histograms.__doc__ = Plots.distributions.__doc__
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 c031aa03f..2bba3a2c2 100644
--- a/OpenPNM/Geometry/__Voronoi__.py
+++ b/OpenPNM/Geometry/__Voronoi__.py
@@ -9,8 +9,16 @@ 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 OpenPNM.Utilities import topology
+logger = logging.getLogger(__name__)
+topo = topology()
 
 
 class Voronoi(GenericGeometry):
@@ -19,15 +27,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,29 +62,297 @@ 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()
 
+        topo.trim_occluded_throats(self._net, self.name)
+
+        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)
         self.models.add(propname='throat.seed',
                         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)
-        self.models.add(propname='throat.volume',
-                        model=gm.throat_volume.extrusion)
+                        const=self._fibre_rad*2)
+        self['throat.volume'] = 0.0
         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]]
+
+        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
+        topo.trim_occluded_throats(self._net)
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..7db1941b5 100644
--- a/OpenPNM/Geometry/models/pore_centroid.py
+++ b/OpenPNM/Geometry/models/pore_centroid.py
@@ -1,37 +1,31 @@
 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]
+    value.fill(0.0)
+    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=True)['target']
         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)]
+        if len(verts) > 0:
+            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..f56e4255a 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=True)['target']
         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..906d63b91 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,29 +345,38 @@ 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')
 
         return A
 
-    def trim_occluded_throats(self):
+    def _export_vor_fibres(self):
         r"""
-        After the offsetting routine throats with zero area have been fully occluded.
-        Remove these from the network and also remove pores that are isolated
+        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
         """
-        occluded_ts = list(self.throats()[self['throat.area'] == 0])
-        if len(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)
+        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 9f9aa5500..37c2c6325 100644
--- a/OpenPNM/Physics/__GenericPhysics__.py
+++ b/OpenPNM/Physics/__GenericPhysics__.py
@@ -101,9 +101,9 @@ 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/capillary_pressure.py b/OpenPNM/Physics/models/capillary_pressure.py
index 077860749..a4043747f 100644
--- a/OpenPNM/Physics/models/capillary_pressure.py
+++ b/OpenPNM/Physics/models/capillary_pressure.py
@@ -41,19 +41,25 @@ def washburn(physics, phase, network, surface_tension='pore.surface_tension',
     suitable for highly non-wetting invading phases in most materials.
 
     """
-    if surface_tension.split('.')[0] == 'pore':
+    if ((surface_tension.split('.')[0] == 'pore') and
+       (throat_diameter.split('.')[0] == 'throat')):
         sigma = phase[surface_tension]
         sigma = phase.interpolate_data(data=sigma)
     else:
         sigma = phase[surface_tension]
-    if contact_angle.split('.')[0] == 'pore':
+    if ((contact_angle.split('.')[0] == 'pore') and
+       (throat_diameter.split('.')[0] == 'throat')):
         theta = phase[contact_angle]
         theta = phase.interpolate_data(data=theta)
     else:
         theta = phase[contact_angle]
     r = network[throat_diameter]/2
     value = -2*sigma*_sp.cos(_sp.radians(theta))/r
-    value = value[phase.throats(physics.name)]
+    if throat_diameter.split('.')[0] == 'throat':
+        value = value[phase.throats(physics.name)]
+    else:
+        value = value[phase.pores(physics.name)]
+    value[_sp.absolute(value) == _sp.inf] = 0
     return value
 
 
@@ -101,12 +107,14 @@ def purcell(physics, phase, network, r_toroid,
     TODO: Triple check the accuracy of this equation
     """
 
-    if surface_tension.split('.')[0] == 'pore':
+    if ((surface_tension.split('.')[0] == 'pore') and
+       (throat_diameter.split('.')[0] == 'throat')):
         sigma = phase[surface_tension]
         sigma = phase.interpolate_data(data=sigma)
     else:
         sigma = phase[surface_tension]
-    if contact_angle.split('.')[0] == 'pore':
+    if ((contact_angle.split('.')[0] == 'pore') and
+       (throat_diameter.split('.')[0] == 'throat')):
         theta = phase[contact_angle]
         theta = phase.interpolate_data(data=theta)
     else:
@@ -117,7 +125,10 @@ def purcell(physics, phase, network, r_toroid,
     value = (-2*sigma/r) * \
         (_sp.cos(_sp.radians(theta - alpha)) /
             (1 + R/r*(1 - _sp.cos(_sp.radians(alpha)))))
-    value = value[phase.throats(physics.name)]
+    if throat_diameter.split('.')[0] == 'throat':
+        value = value[phase.throats(physics.name)]
+    else:
+        value = value[phase.pores(physics.name)]
     return value
 
 
diff --git a/OpenPNM/Physics/models/diffusive_conductance.py b/OpenPNM/Physics/models/diffusive_conductance.py
index b88fab87e..e8674e29a 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,15 @@ 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())
+    sf[_sp.isnan(sf)] = 1.0
+    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..f4a16c12f 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,13 @@ 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())
+    sf[_sp.isnan(sf)] = 1.0
+    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 377ebc9f3..50c9d6ad2 100644
--- a/OpenPNM/Physics/models/multiphase.py
+++ b/OpenPNM/Physics/models/multiphase.py
@@ -100,7 +100,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 0fb6b20f7..54bd9204f 100644
--- a/OpenPNM/Postprocessing/Plots.py
+++ b/OpenPNM/Postprocessing/Plots.py
@@ -178,35 +178,36 @@ def saturation_profile(network, phase, fig=None, axis=2):
     return fig
 
 
-def distributions(net, fig=None, throat_diameter='throat.diameter',
-                  pore_diameter='pore.diameter', throat_length='throat.length',
-                  exclude_boundaries=True, geom_list=None):
+def distributions(obj,
+                  throat_diameter='throat.diameter',
+                  pore_diameter='pore.diameter',
+                  throat_length='throat.length'):
     r"""
     Plot a montage of key network size distribution histograms
 
     Parameters
     ----------
-    net : OpenPNM Network Object
-    The network for which the graphs are desired
+    obj : OpenPNM Object
+        This object can either be a Network or a Geometry.  If a Network is
+        sent, then the histograms will display the properties for the entire
+        Network.  If a Geometry is sent, then only it's properties will be
+        shown.
+    throat_diameter : string
+        Dictionary key to the array containing throat diameter values
+    pore_diameter : string
+        Dictionary key to the array containing pore diameter values
+    throat_length : string
+        Dictionary key to the array containing throat length values
 
     """
-    if fig is None:
-        fig = _plt.figure()
+    fig = _plt.figure()
 
     fig.subplots_adjust(hspace=0.4)
     fig.subplots_adjust(wspace=0.4)
 
-    if geom_list is not None:
-        include_pores = [False]*net.num_pores()
-        include_throats = [False]*net.num_throats()
-        for geom in geom_list:
-            include_pores = include_pores | net['pore.' + geom]
-            include_throats = include_throats | net['throat.' + geom]
-    else:
-        include_pores = net['pore.all']
-        include_throats = net['throat.all']
-    pores = net.pores()[include_pores]
-    throats = net.throats()[include_throats]
+    pores = obj._net.pores(obj.name)
+    throats = obj._net.throats(obj.name)
+    net = obj._net
 
     ax1 = fig.add_subplot(221)
     ax1.hist(net[pore_diameter][pores], 25, facecolor='green')
diff --git a/OpenPNM/Utilities/__topology__.py b/OpenPNM/Utilities/__topology__.py
index ec5aeafed..ea8a17152 100644
--- a/OpenPNM/Utilities/__topology__.py
+++ b/OpenPNM/Utilities/__topology__.py
@@ -446,7 +446,7 @@ class topology(object):
 
         Parameters
         ----------
-        networK : OpenPNM Network Object
+        network : OpenPNM Network Object
 
         pores : array_like
             The first group of pores to be replaced
@@ -583,3 +583,29 @@ class topology(object):
         for l in main_labels:
             del network['pore.surface_'+l]
         self.trim(network=network, pores=pores)
+
+    def trim_occluded_throats(self, network, mask='all'):
+        r"""
+        Remove throats with zero area from the network and also remove
+        pores that are isolated (as a result or otherwise)
+
+        Parameters
+        ----------
+        network : OpenPNM Network Object
+
+        mask : string
+            Applies routine only to pores and throats with this label
+        """
+        occluded_ts = network['throat.area'] == 0
+        if _sp.sum(occluded_ts) > 0:
+            # Apply mask
+            occluded_ts *= network["throat."+mask]
+            self.trim(network=network, throats=occluded_ts)
+        # Also get rid of isolated pores
+        isolated_ps = network.check_network_health()['isolated_pores']
+        if _sp.size(isolated_ps) > 0:
+            # Convert to Bool array and apply mask
+            temp_array = _sp.zeros(network.num_pores()).astype(bool)
+            temp_array[isolated_ps] = True
+            isolated_ps = temp_array * network["pore."+mask]
+            self.trim(network=network, pores=isolated_ps)
diff --git a/OpenPNM/Utilities/misc.py b/OpenPNM/Utilities/misc.py
index c2a0188e5..ac32c116a 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..65dff4e0d 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 plot_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-')
-            fig.show()
+            ax.ticklabel_format(style='sci', scilimits=(0, 0))
+
     else:
-        print('Please provide throat indices')
+        print("Please provide throat indices")
+    return fig
 
 
-def print_pore(geom, pores, fig=None, axis_bounds=None):
+def plot_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')
-            plt.show()
+            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))
+
         else:
-            print_throat(throats)
+            plot_throat(geometry, throats, fig)
     else:
         print('Please provide pore indices')
-    if return_fig is True:
-        return fig
+    return fig
diff --git a/docs/examples/voronoi.rst b/docs/examples/voronoi.rst
index 70d5cc220..758d69087 100644
--- a/docs/examples/voronoi.rst
+++ b/docs/examples/voronoi.rst
@@ -10,7 +10,7 @@ Random PNM with Delaunay Network and Voronoi Geometry
 Getting Started
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-In this tutorial, we will demonstrate concepts of random pore network modelling outlined in J.T. Gostick's 2013 paper `[1]`_. Pores are randomly dispersed inside the domain and connections are defined by the nearest neighbour Delaunay triangulation algorithm. The Voronoi diagram is the copliment to the Delaunay trinagulation and is used to replicate the fibres in a porous domain such as a fuel cell GDL. Each pore is surrounded by a cage of fibres which are created from the intersections of the equidistant planes between neighbouring pores. Where the planes intersect, a Voronoi vertex is created, and these are saved to the Delaunay Network object as "pore.vertices". A throat is defined by a set of vertices which are shared by two neighbouring pores. The throat vertices will be coplanar and the throat normal is the vector from one pore's co-ordinate to the other. The vertices are used by the Voronoi geometry model which effectively erodes the fibre space by offsetting the vertices in the planes of each connecting throat by a fibre radius parameter passed to the method.
+In this tutorial, we will demonstrate concepts of random pore network modelling outlined in J.T. Gostick's 2013 paper `[1]`_. Pores are randomly dispersed inside the domain and connections are defined by the nearest neighbour Delaunay triangulation algorithm. The Voronoi diagram is the compliment to the Delaunay trinagulation and is used to replicate the fibres in a porous domain such as a fuel cell GDL. Each pore is surrounded by a cage of fibres which are created from the intersections of the equidistant planes between neighbouring pores. Where the planes intersect, a Voronoi vertex is created, and these are saved to the Delaunay Network object as "pore.vert_index" and "throat.vert_index". A throat is defined by a set of vertices which are shared by two neighbouring pores. The throat vertices will be coplanar and the throat normal is the vector from one pore's co-ordinate to the other. The vertices are used by the Voronoi geometry model which creates a 3D image of the fibres using a supplied fibre radius. Image analysis is then performed to extract pore and throat sizes by using the convex hull of the pore and throat vertices.
  
 .. _[1]: http://jes.ecsdl.org/cgi/doi/10.1149/2.009308jes
 
@@ -22,52 +22,51 @@ There will be a general layout to complete this simulation:
 4. Set up phase physics and methods 
 5. Run ordinary percolation 
 6. Plot a drainage curve
+   
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Setting up Network and Geometry
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 We first import the OpenPNM code.
 
 .. code-block:: python
     
     import OpenPNM
-   
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Setting up Network and Geometry
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 To begin our simulation, we must first generate our random network and geometry.  This includes:
 
-1. Creating a Delaunay network object and a Voronoi geometry object
-2. Creating boundary pores
+1. Getting a handle to the OpenPNM Controller and setting the loglevel to display useful information
+2. Creating a Delaunay network object and a Voronoi geometry object
 3. Assingning a 'Voronoi' geometry to the initial set of pores and a 'Boundary' geometry to the boundary pores
 
 .. code-block:: python
 
-	pn = OpenPNM.Network.Delaunay(num_pores=100, domain_size=[1,1,1],name='net')
+	ctrl = OpenPNM.Base.Controller()
+	ctrl.loglevel = 20
+	pn = OpenPNM.Network.Delaunay(num_pores=100, domain_size=[1e-4, 1e-4, 1e-4], name='net')
 	pn.add_boundaries()
 	# create geometry
-	fibre_rad=0.03
-	Ps = pn.pores('boundary',mode='not')
-	Ts = pn.find_neighbor_throats(pores=Ps,mode='intersection',flatten=True)
-	geo = OpenPNM.Geometry.Voronoi(network=pn,pores=Ps,throats=Ts,fibre_rad=fibre_rad)
-	# create boundary geometry
-	Ps = pn.pores('boundary')
-	Ts = pn.find_neighbor_throats(pores=Ps,mode='not_intersection')
-	boun = OpenPNM.Geometry.Boundary(network=pn,pores=Ps,throats=Ts,name='boun')
+	fibre_rad = 5e-6
+	Ps = pn.pores()
+	Ts = pn.find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)
+	geo = OpenPNM.Geometry.Voronoi(network=pn, pores=Ps, throats=Ts, fibre_rad=fibre_rad, voxel_vol=True, vox_len=1e-6, name='vor')
 	
-It is worth mentioning a little about the boundaries at this point. Behind the scenes all the pores in the network were temporarily reflected about the planes confining the pore space during the Delaunay tesselation. This has the effect of creating throats on the outer confines that are aligned with the boundary planes. The add_boundaries() method makes use of this feature and generates 6 new sets of boundary pores which have their co-ordinates placed on the boundary plane in-line with the adjoining pore. The boundary pores are labelled "left_boundary", "right_boundary", etc. and all boundary pores are also given the label "boundary". It is generally a good idea to assign these pores to a Boundary type Geometry which will give the pores zero volume.
+It is worth mentioning a little about the boundaries at this point. Behind the scenes all the pores in the network were temporarily reflected about the planes confining the pore space during the Delaunay tesselation. This has the effect of creating throats on the outer confines that are aligned with the boundary planes. The add_boundaries() method makes use of this feature and generates 6 new sets of boundary pores which have their co-ordinates placed on the boundary plane in-line with the adjoining pore. The boundary pores are labelled "left_boundary", "right_boundary", etc. and all boundary pores are also given the label "boundary".
 	
-The 'Voronoi' geometry object has two visualisation methods which are unique to it. Both methods will plot a selection of pores or throats and help visualise and check the offsetting routine.  
+The 'vertexops' utility has two visualisation methods which are unique to Voronoi geometries. Both methods will plot a selection of pores or throats and help visualise and check the offsetting routine.  
 
 .. code-block:: python
 
-	#print all pores in the Voronoi geometry
-	geo.print_pore(geo.map_pores(target=pn,pores=geo.Ps))
+	#plot all pores in the Voronoi geometry
+	import OpenPNM.Utilities.vertexops as vo
+	vo.plot_pore(geo, geo.pores())
 	throats = pn.find_neighbor_throats(pores=[0])
-	#print all throats connected to the first pore in the network
-	geo.print_throat(throats)
+	#plot all throats connected to the first pore in the network
+	vo.plot_throat(geo, throats)
 
-.. image:: http://imgur.com/zrEFZ1O.png
+.. image:: http://imgur.com/icZfxQ2.png
 
-.. image:: http://imgur.com/mmnQsHz.png
+.. image:: http://imgur.com/sCL1kY1.png
 
 We can also see the distributions of the pore and throat sizes.
 
@@ -75,9 +74,91 @@ We can also see the distributions of the pore and throat sizes.
 
 	OpenPNM.Postprocessing.Plots.distributions(net=pn)
 
-.. image:: http://imgur.com/Zw93uSV.png
+.. image:: http://imgur.com/IylwVk5.png
+
+The Voronoi geometry has two parameters which control whether a voxel image of the fibres is generated or not and what resolution is used. If 'voxel_vol' is set to True then the image is created with resolution set by the 'vox_len' parameter. The default is to generate the image at resolution of 1e-6 voxel side length giving each voxel a volume of 1e-18. Generating the voxel image is a memory intensive process relying on many image analysis routines and it is recommended that a smaller network is tested first on your machine whilst monitroing your system performance to guage whether larger networks are possible. Setting 'voxel_vol' to False will mean that pore volumes are calculated from the ConvexHull of the offset throat vertices which is a faster but less accurate method. 
+
+If the voxel image is generated a few more methods of the Voronoi Geometry class can be used to visualize and export the image.
+
+.. code-block:: python
+
+	geo.plot_fibre_slice(plane=[0.5,0,0])
+	
+.. image:: http://imgur.com/GMzeBck.png
+
+.. code-block:: python
+
+	geo.plot_porosity_profile()
+
+.. image:: http://imgur.com/j0kTCRj.png
 
-Three new pieces of information are added to the Network object when the Vornoi geometry is applied to a set of pores: 'throat.verts', 'throat.offset_verts' and 'throat.normals'. The offset vertices are subsequently used to calculate the throat area. Also the pore volume is calculated using the convex hull of all the offset throat vertices for a given pore.
+The following command creates a binary array which can be used in Matlab.
+
+.. code-block:: python
+
+	geo._export_fibre_image()
+
+In addition to visualization within OpenPNM using the voxel image it is possible to export a pickle dump of the throat vertices in convex hull order which form the skeleton of the fibres. 
+
+.. code-block:: python
+
+	pn._export_vor_fibres()
+
+This can then be used in Blender to create images such as this:
+
+.. image:: http://imgur.com/MCPQ1S7.png
+
+The Blender script used to create the multi-faceted Tetrahedron is as follows:
+
+.. code-block:: python
+
+	import bpy
+	import mathutils
+	import numpy as np
+	import pickle as pickle
+	 
+	class TetrahedronMakerPanel(bpy.types.Panel):
+		bl_space_type = "VIEW_3D"
+		bl_region_type = "TOOLS"
+		bl_context = "objectmode"
+		bl_category = "Create"
+		bl_label = "Add Tetrahedron"
+	 
+		def draw(self, context):
+			TheCol = self.layout.column(align=True)
+			TheCol.operator("mesh.make_tetrahedron", text="Add Tetrahedron")
+		#end draw
+	 
+	#end TetrahedronMakerPanel
+	 
+	class MakeTetrahedron(bpy.types.Operator):
+		bl_idname = "mesh.make_tetrahedron"
+		bl_label = "Add Tetrahedron"
+		bl_options = {"UNDO"}
+	 
+		def invoke(self, context, event):
+			data = pickle.load( open( "C:\\Working_Files\\OpenPNM\\fibres.p", "rb" ) )
+			Verts = np.around(data["Verts"]*1e4,2)
+			Faces = data['Indices']
+			NewMesh = bpy.data.meshes.new("Tetrahedron")
+			NewMesh.from_pydata \
+			  (
+				Verts,
+				[],
+				Faces
+			  )
+			NewMesh.update()
+			NewObj = bpy.data.objects.new("Tetrahedron", NewMesh)
+			context.scene.objects.link(NewObj)
+			return {"FINISHED"}
+		#end invoke
+	 
+	#end MakeTetrahedron
+	 
+	bpy.utils.register_class(MakeTetrahedron)
+	bpy.utils.register_class(TetrahedronMakerPanel)
+
+The image is then created by deleteing the faces leaving just the wire-frame edges of the Voronoi diagram and then converting the object to a curve and applying a circular bevel curve.
 
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Setting up the Phases and Physics
@@ -91,10 +172,7 @@ Now we are ready to set up our phases (water and air) and the physics correspond
     air = OpenPNM.Phases.Air(network = pn, name = 'air')
     water = OpenPNM.Phases.Water(network = pn, name = 'water')
 
-We are now ready to establish physical properties for our fluid objects. To do this, we will: 
-
-1. Create physics objects associated with our fluids (by using BasePhyics we don't have to add methods for calculating each property because they are already included) 
-2. Use our regenerate_physics() method to calculate these properties
+We are now ready to establish physical properties for our fluid objects. To do this, we will create physics objects associated with our fluids (by using OpenPNM.Physics.Standard we don't have to add methods for calculating each property because they are already included) 
 
 .. code-block:: python
 
@@ -114,7 +192,7 @@ To run our simulation, we will:
 1. Pick inlet pores
 2. Create an Ordinary Percolation algorithm object
 3. Run our algorithm object
-4. Call update() so that occupancy of pores and throats for each fluid will be set
+4. Return results so that occupancy of pores and throats for each fluid will be set
 
 .. code-block:: python
 
@@ -125,8 +203,8 @@ To run our simulation, we will:
     #prevents extremely small diffusivity and permeability values in the z direction
     used_inlets = [inlets[x] for x in range(0, len(inlets), 2)]
     
-    OP_1 = OpenPNM.Algorithms.OrdinaryPercolation(network=pn,invading_phase=water,defending_phase=air)
-    OP_1.run(inlets=used_inlets,npts=100)
+    OP_1 = OpenPNM.Algorithms.OrdinaryPercolation(network=pn, invading_phase=water, defending_phase=air)
+    OP_1.run(inlets=used_inlets, npts=100)
 
 This algorithm performed a start to finish simulation, which fully flooded the network. The 'return_results()' command can be used to update the phase occupancy values throughout the network. 
 
@@ -139,7 +217,7 @@ OpenPNM makes it very easy to inspect the ouput of the algorithm through the "Po
 
 .. code-block:: python
 
-	OpenPNM.Postprocessing.Plots.drainage_curves(OP_1,timing=None)
+	OpenPNM.Postprocessing.Plots.drainage_curves(OP_1, timing=None)
 
 We can also view the network data by creating vtk files to be viewed using ParaView (downloadable at http://www.paraview.org/download/ ). It is suggested that version 3.98 is downloaded instead of 4.1).  If we visualize our pore network model with phase data included it will look like this:
 
@@ -151,14 +229,7 @@ To create the vtk file use the following command
 
 .. code-block:: python
 	
-    import OpenPNM.Utilities.IO as io
-    io.VTK.save(network=pn,phases=[air,water])
-
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Differences between OpenPNM and Gostick's simulation
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-The major difference between the technique described in Gostick's paper `[1]`_ and that implemented in OpenPNM is the method used to replicate the fibres. Where Gostick uses morphological image analysis to construct 3D voxel images of the solid fibres a simpler method using geometrical offsetting of the Voronoi vertices is used in OpenPNM. This technique is quicker than image analysis and can be achieved entirely within the framework of OpenPNM, however, a fibrous pore space is less well represented as fibres are effectively treated as polygons rather than cylinders.
+    ctrl.export(pn)
 
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 References
diff --git a/run_script.py b/run_script.py
index 2ae6c4308..e404c2c74 100644
--- a/run_script.py
+++ b/run_script.py
@@ -1,6 +1,7 @@
 import OpenPNM as op
-print('-----> Using OpenPNM version: '+op.__version__)
-
+import time
+st = time.time()
+from OpenPNM.Geometry import models as gm
 #==============================================================================
 '''Build Topological Network'''
 #==============================================================================
@@ -13,7 +14,7 @@ pn.add_boundaries()
 Ps = pn.pores('boundary',mode='not')
 Ts = pn.find_neighbor_throats(pores=Ps,mode='intersection',flatten=True)
 geom = op.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 = op.Geometry.Boundary(network=pn,pores=Ps,throats=Ts)
@@ -92,3 +93,5 @@ except Exception as e:
 '''Export to VTK'''
 #------------------------------------------------------------------------------
 op.export()
+
+print("sim time:" + str(time.time()-st))
\ No newline at end of file
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..16a4d640c 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..be3dee578 100644
--- a/test/unit/Network/DelaunayTest.py
+++ b/test/unit/Network/DelaunayTest.py
@@ -1,3 +1,31 @@
+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_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/Postprocessing/PlotsTest.py b/test/unit/Postprocessing/PlotsTest.py
index d84390e72..1ededd1c5 100644
--- a/test/unit/Postprocessing/PlotsTest.py
+++ b/test/unit/Postprocessing/PlotsTest.py
@@ -4,6 +4,8 @@ import matplotlib as mpl
 
 class PlotsTest:
     def setup_class(self):
+        ctrl = OpenPNM.Base.Controller()
+        ctrl.clear()
         self.net = OpenPNM.Network.Cubic(shape=[10, 10, 10])
         self.geo = OpenPNM.Geometry.Toray090(network=self.net,
                                              pores=self.net.Ps,
diff --git a/test/unit/Utilities/TopologyTest.py b/test/unit/Utilities/TopologyTest.py
index 906829faf..b62ccc046 100644
--- a/test/unit/Utilities/TopologyTest.py
+++ b/test/unit/Utilities/TopologyTest.py
@@ -1,3 +1,28 @@
+import OpenPNM
+from OpenPNM.Utilities import topology
+import scipy as sp
+topo = topology()
+
+
 class TopologyTest:
+
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5], spacing=1)
+        Ps = self.net.pores()
+        Ts = self.net.throats()
+        self.geom = OpenPNM.Geometry.TestGeometry(network=self.net,
+                                                  pores=Ps,
+                                                  throats=Ts,
+                                                  name='test_geom')
+
     def test_extend(self):
         pass
+
+    def test_trim_occluded_throats(self):
+        Np = self.net.num_pores()
+        Nt = self.net.num_throats()
+        ts = self.net.find_neighbor_throats(pores=[0])
+        self.geom['throat.area'][ts] = 0.0
+        topo.trim_occluded_throats(self.net)
+        assert self.net.num_pores() == Np - 1
+        assert self.net.num_throats() == Nt - sp.size(ts)
diff --git a/test/unit/Utilities/VertexOpsTest.py b/test/unit/Utilities/VertexOpsTest.py
index 73ae451da..02166e5bc 100644
--- a/test/unit/Utilities/VertexOpsTest.py
+++ b/test/unit/Utilities/VertexOpsTest.py
@@ -1,3 +1,56 @@
+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
+
+    def test_plot_pore(self):
+        vo.plot_pore(self.geo, self.geo.pores())
+
+    def test_plot_throat(self):
+        vo.plot_throat(self.geo, [1, 2, 3, 4])