diff --git a/openpnm/_skgraph/generators/_voronoi_delaunay_dual.py b/openpnm/_skgraph/generators/_voronoi_delaunay_dual.py index fac3c1fd1b..3ddd9cf070 100644 --- a/openpnm/_skgraph/generators/_voronoi_delaunay_dual.py +++ b/openpnm/_skgraph/generators/_voronoi_delaunay_dual.py @@ -22,11 +22,14 @@ def voronoi_delaunay_dual(points, shape, trim=True, reflect=True, trim : bool, optional If ``True`` (default) then all points lying beyond the given domain shape will be removed + reflect : bool, optionl + If ``True`` (Default) then points are reflected across each face of the + domain prior to performing the tessellation. Returns ------- network : dict - A dictionary containing 'node.coords' and 'edge.conns' + A dictionary containing '.coords' and '.conns' vor : Voronoi object The Voronoi tessellation object produced by ``scipy.spatial.Voronoi`` tri : Delaunay object @@ -110,18 +113,22 @@ def voronoi_delaunay_dual(points, shape, trim=True, reflect=True, Ts = np.sum(network[node_prefix+'.delaunay'][conns].astype(int), axis=1) == 1 network[edge_prefix+'.interconnect'] = Ts - # Identify and trim nodes outside the domain if requested if trim: - inside_all = ~isoutside(network, shape=shape) - inside_delaunay = inside_all*network[node_prefix+'.delaunay'] - outside_delaunay = (~inside_all)*network[node_prefix+'.delaunay'] - neighbors = find_neighbor_nodes(network=network, - inds=np.where(inside_delaunay)[0], - include_input=True) - trim = np.ones([network[node_prefix+'.coords'].shape[0], ], dtype=bool) - trim[neighbors] = False # Keep all neighbors to internal delaunay nodes - trim[outside_delaunay] = True # Re-add external delaunay nodes to trim - network = trim_nodes(network=network, inds=np.where(trim)[0]) + # Find all delaunay nodes outside the domain + Ps = isoutside(network=network, shape=shape)*network[node_prefix+'.delaunay'] + if np.any(Ps): # only occurs if points were reflected + # Find voronoi nodes connected to these and mark them as surface nodes + inds = np.where(Ps)[0] + Ns = find_neighbor_nodes(network=network, inds=inds) + network[node_prefix+'.surface'] = np.zeros(n_nodes, dtype=bool) + network[node_prefix+'.surface'][Ns] = True + Ps = isoutside(network=network, shape=shape) + inds = np.where(Ps)[0] + network = trim_nodes(network=network, inds=inds) + else: + trim = isoutside(network=network, shape=shape) + inds = np.where(trim)[0] + network = trim_nodes(network=network, inds=inds) return network, vor, tri diff --git a/openpnm/network/_delaunay.py b/openpnm/network/_delaunay.py index 5dfc34a421..5b95ddb31b 100644 --- a/openpnm/network/_delaunay.py +++ b/openpnm/network/_delaunay.py @@ -1,6 +1,6 @@ from openpnm.network import Network from openpnm.utils import Docorator -from openpnm._skgraph.generators import delaunay, tools +from openpnm._skgraph.generators import delaunay __all__ = ['Delaunay'] @@ -39,12 +39,6 @@ class Delaunay(Network): %(Network.parameters)s - See Also - -------- - Gabriel - Voronoi - DelaunayVoronoiDual - Notes ----- It is also possible to generate circular ``[r, 0]``, cylindrical ``[r, z]``, and @@ -59,7 +53,6 @@ class Delaunay(Network): def __init__(self, shape, points, reflect=True, trim=True, **kwargs): super().__init__(**kwargs) - points = tools.parse_points(shape=shape, points=points, reflect=reflect) net, tri = delaunay(points=points, shape=shape, node_prefix='pore', diff --git a/openpnm/network/_delaunay_voronoi_dual.py b/openpnm/network/_delaunay_voronoi_dual.py index 93d8ce3448..75b96e5fc4 100644 --- a/openpnm/network/_delaunay_voronoi_dual.py +++ b/openpnm/network/_delaunay_voronoi_dual.py @@ -1,14 +1,15 @@ import numpy as np -import scipy.spatial as sptl -from openpnm import topotools from openpnm.network import Network from openpnm._skgraph.generators.tools import parse_points from openpnm._skgraph.generators import voronoi_delaunay_dual +from openpnm.utils import Docorator __all__ = ['DelaunayVoronoiDual'] +docstr = Docorator() +@docstr.dedent class DelaunayVoronoiDual(Network): r""" Combined and interconnected Voronoi and Delaunay tessellations @@ -19,55 +20,48 @@ class DelaunayVoronoiDual(Network): Can either be an N-by-3 array of point coordinates which will be used, or a scalar value indicating the number of points to generate shape : array_like - The size and shape of the domain used for generating and trimming - excess points. The coordinates are treated as the outer corner of a - rectangle [x, y, z] whose opposite corner lies at [0, 0, 0]. - By default, a domain size of [1, 1, 1] is used. To create a 2D network - set the Z-dimension to 0. + The size and shape of the domain: + + ========== ============================================================ + shape result + ========== ============================================================ + [x, y, z] A 3D cubic domain of dimension x, y and z + [x, y, 0] A 2D square domain of size x by y + ========== ============================================================ + + trim : bool, optional + If ``True`` (default) then all vertices laying outside the domain will + be removed. This is only useful if ``reflect=True``. + reflect : bool, optional + If ``True`` (default) then the base points will be reflected across + all the faces of the domain prior to performing the tessellation. This + feature is best combined with ``trim=True`` to prevent unreasonably long + connections between points on the surfaces. + + %(Network.parameters)s Notes ----- A Delaunay tessellation is performed on a set of base points then the corresponding Voronoi diagram is generated. Finally, each Delaunay node - is connected to it's neighboring Voronoi vertices to create interaction + is connected to its neighboring Voronoi vertices to create interconnections between the two networks. All pores and throats are labelled according to their network (i.e. 'pore.delaunay'), so they can be each assigned to a different Geometry. - The dual-nature of this network is meant for modeling transport in the void - and solid space simultaneously by treating one network (i.e. Delaunay) as - voids and the other (i.e. Voronoi) as solid. Interaction such as heat - transfer between the solid and void can be accomplished via the - interconnections between the Delaunay and Voronoi nodes. - """ def __init__(self, shape, points, trim=True, reflect=True, **kwargs): super().__init__(**kwargs) - points = parse_points(shape=shape, points=points, reflect=reflect) - net, vor, tri = voronoi_delaunay_dual(shape=shape, points=points, trim=trim, node_prefix='pore', edge_prefix='throat') self.update(net) - self._vor = vor - self._tri = tri - - @property - def tri(self): - """A shortcut to get a handle to the Delanuay subnetwork""" - if not hasattr(self, '_tri'): - points = self._vor.points - self._tri = sptl.Delaunay(points=points) - return self._tri - - @property - def vor(self): - """A shortcut to get a handle to the Voronoi subnetwork""" - return self._vor + self.vor = vor + self.tri = tri def find_throat_facets(self, throats=None): r""" diff --git a/openpnm/network/_voronoi.py b/openpnm/network/_voronoi.py index 5ddaad6e74..cc9001a68b 100644 --- a/openpnm/network/_voronoi.py +++ b/openpnm/network/_voronoi.py @@ -1,7 +1,6 @@ from openpnm.network import Network from openpnm.utils import Docorator from openpnm._skgraph.generators import voronoi -from openpnm._skgraph.generators.tools import parse_points docstr = Docorator() @@ -60,14 +59,10 @@ class Voronoi(Network): def __init__(self, shape, points, trim=True, reflect=True, **kwargs): super().__init__(**kwargs) - # Clean-up input points - points = parse_points(shape=shape, points=points, reflect=reflect) - # Call skgraph generator net, vor = voronoi(points=points, shape=shape, trim=trim, node_prefix='pore', edge_prefix='throat') - # Update dictionary on self self.update(net) self.vor = vor diff --git a/tests/unit/network/VoronoiDelaunayDualTest.py b/tests/unit/network/VoronoiDelaunayDualTest.py new file mode 100644 index 0000000000..0fe404e18d --- /dev/null +++ b/tests/unit/network/VoronoiDelaunayDualTest.py @@ -0,0 +1,107 @@ +import numpy as np +import openpnm as op + + +class VoronoiTest: + + def setup_class(self): + pass + + def teardown_class(self): + pass + + def test_dual_square_trim_and_reflect(self): + np.random.seed(0) + shape = [1, 1, 0] + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=False, reflect=False) + assert op.topotools.isoutside(network=net, shape=shape).sum() > 0 + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=True, reflect=False) + assert op.topotools.isoutside(network=net, shape=shape).sum() == 0 + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=False, reflect=True) + assert op.topotools.isoutside(network=net, shape=shape).sum() > 0 + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=True, reflect=True) + assert op.topotools.isoutside(network=net, shape=shape).sum() == 0 + + def test_dual_cube_trim_and_reflect(self): + np.random.seed(0) + shape = [1, 1, 1] + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=False, reflect=False) + assert op.topotools.isoutside(network=net, shape=shape).sum() > 0 + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=True, reflect=False) + assert op.topotools.isoutside(network=net, shape=shape).sum() == 0 + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=False, reflect=True) + assert op.topotools.isoutside(network=net, shape=shape).sum() > 0 + net = op.network.DelaunayVoronoiDual(points=30, shape=shape, + trim=True, reflect=True) + assert op.topotools.isoutside(network=net, shape=shape).sum() == 0 + + def test_dual_square_num_points(self): + np.random.seed(0) + net = op.network.DelaunayVoronoiDual(points=30, shape=[1, 1, 0]) + assert net.Np > 30 + assert np.all(net.coords[:, -1] == 0.0) + assert np.all(op.topotools.dimensionality(net) == [True, True, False]) + + def test_dual_cube_num_points(self): + np.random.seed(0) + net = op.network.DelaunayVoronoiDual(points=30, shape=[1, 1, 1]) + assert net.Np > 30 + assert not np.all(net.coords[:, -1] == 0.0) + assert np.all(op.topotools.dimensionality(net) == [True, True, True]) + + def test_dual_disk_num_points(self): + np.random.seed(0) + net = op.network.DelaunayVoronoiDual(points=30, shape=[1, 0]) + assert np.all(net.coords[:, -1] == 0.0) + + def test_dual_cylinder_num_points(self): + np.random.seed(0) + net = op.network.DelaunayVoronoiDual(points=30, shape=[1, 1]) + assert not np.all(net.coords[:, -1] == 0.0) + assert np.all(op.topotools.dimensionality(net) == [True, True, True]) + + def test_dual_sphere_num_points(self): + np.random.seed(0) + net = op.network.DelaunayVoronoiDual(points=30, shape=[1]) + assert not np.all(net.coords[:, -1] == 0.0) + assert np.all(op.topotools.dimensionality(net) == [True, True, True]) + + def test_dual_delaunay_dual_check_points(self): + np.random.seed(0) + dual = op.network.DelaunayVoronoiDual(points=100, shape=[1, 1, 1]) + assert dual.num_pores('delaunay') == 100 + assert dual.num_pores('voronoi') == 567 + + def test_find_throat_facets(self): + np.random.seed(0) + dual = op.network.DelaunayVoronoiDual(points=10, shape=[1, 1, 1]) + f = dual.find_throat_facets(throats=[1, 5]) + assert np.all(f[0] == [48, 49, 50, 55, 57]) + assert np.all(f[1] == [48, 33, 30, 49]) + + def test_find_pore_hulls(self): + np.random.seed(0) + dual = op.network.DelaunayVoronoiDual(points=10, shape=[1, 1, 1]) + f = dual.find_pore_hulls(pores=[0, 5]) + assert np.all(f[0] == [12, 14, 15, 19, 20, 21, 30, + 33, 35, 48, 49, 50, 55, 57]) + assert np.all(f[1] == [36, 37, 38, 39, 40, 41, 42, + 43, 51, 58, 60, 61]) + + +if __name__ == '__main__': + + t = VoronoiTest() + t.setup_class() + self = t + for item in t.__dir__(): + if item.startswith('test'): + print(f'Running test: {item}') + t.__getattribute__(item)()