Adding boundary pores#
In this example we will illustrate how boundary pores are added during the network extraction phase. We will use the individual steps of the SNOW algorithm (i.e.
porespy.filters.snow_partioning), rather than use the pre-packaged
import porespy as ps import matplotlib.pyplot as plt import numpy as np import openpnm as op from edt import edt from copy import copy ws = op.Workspace() ws.settings['loglevel'] = 50 ps.visualization.set_mpl_style()
/opt/hostedtoolcache/Python/3.8.16/x64/lib/python3.8/site-packages/openpnm/algorithms/_invasion_percolation.py:358: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details. def _find_trapped_pores(inv_seq, indices, indptr, outlets): # pragma: no cover
im = ps.generators.lattice_spheres(shape=[52, 52], r=10, offset=13, spacing=25) ps.imshow(im);
snow = ps.filters.snow_partitioning(im) fig, ax = plt.subplots(1, 3, figsize=[10, 4]) ax.imshow(snow.dt, origin='lower', interpolation='none') ax.imshow(snow.peaks > 0, origin='lower', interpolation='none') ax.imshow(snow.regions, origin='lower', interpolation='none');
bds = ps.networks.add_boundary_regions(regions=snow.regions, pad_width=[5, 5]) cm = copy(plt.cm.nipy_spectral) plt.figure(figsize=[4, 4]) plt.imshow(bds, origin='lower', interpolation='none', cmap=cm);
The image above shows how boundary regions appear. They are 3-voxel thick slices that are padded onto each face and given a unique label to differentiate them from the internal pores to which they connect. There are a few key things to note:
Boundary regions that touch each other are separated by a small solid obstacle indicated by black. This is to prevent boundary pores from becoming connected to each other when the network extraction function is applied (
Corners are also converted to solid obstacles. This is for the same reason as above, but also because it’s not clear which internal pores these corners should be connected to.
A thickness of 3 voxels was chosen since this was necessary for the
marching_cubessurface area calculation. It is also a good idea for the pore centroid to be padded by 1-voxel on each side. This 3-voxel dimension can pose problems during the network extraction stage since most of the physical dimensions will be meaningless. Consequently care must be taken to ensure that these boundary pores are treated correctly. In many cases it may be preferrable to forego the them entirely and just apply boundary conditions on real pores lying on the domain surface.
Passing this image into the
regions_to_network function will treat each region as an actual pore:
net = ps.networks.regions_to_network(bds)
Next we’ll open the
net dictionary in OpenPNM using the
PoreSpy IO class:
pn = op.io.network_from_porespy(net) fig = plt.figure() fig = op.visualization.plot_connections(pn, c='w', linewidth=2, ax=fig) fig = op.visualization.plot_coordinates(pn, c='w', s=200, ax=fig) plt.imshow(bds.T, cmap=cm, origin='lower') plt.axis('off');
The approach outlined above is missing one important feature. The boundary pores are not labelled. When the process is performed internally by the
regions_to_network function this is taken care of.