fixes and cleanups of voronoi/delaunay generators, added lots of tests

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 '<node_prefix>.coords' and '<edge_prefix>.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
 

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',

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"""

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

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)()