New release (minor): v3.3.0

.github/workflows/examples.yml

@@ -36,6 +36,7 @@ jobs:
           pip install \
             -r requirements.txt \
             -r requirements/tests.txt
+          pip install pypardiso
 
       - name: Running tests
         run:

.github/workflows/gh-pages.yml

@@ -36,6 +36,7 @@ jobs:
         run: |
           pip install -r requirements.txt
           pip install -r requirements/docs.txt
+          pip install pypardiso
 
       - name: Build the documentation
         run: |

.github/workflows/nightly.yml

@@ -49,6 +49,11 @@ jobs:
             -r requirements.txt \
             -r requirements/tests.txt
 
+      - name: Install pypardiso on non-macOS
+        if: (matrix.os != 'macos-latest')
+        run: |
+          pip install pypardiso            
+
       - name: Running tests
         run:
           pytest . --color=yes

.github/workflows/tests.yml

@@ -46,6 +46,11 @@ jobs:
             -r requirements.txt \
             -r requirements/tests.txt
 
+      - name: Install pypardiso on non-macOS
+        if: (matrix.os != 'macos-latest')
+        run: |
+          pip install pypardiso
+
       - name: Disable numba JIT for codecov to include jitted methods
         if: (matrix.python-version == 3.9) && (matrix.os == 'ubuntu-latest')
         run: |

openpnm/__init__.py

@@ -8,25 +8,34 @@
 
 """
 
-from . import _skgraph
-from . import utils
-from . import core
-from . import models
-from . import topotools
-from . import network
-from . import phase
-from . import algorithms
-from . import solvers
-from . import integrators
-from . import io
-from . import contrib
-from . import visualization
-
-from .utils import Workspace, Project
+import logging
+
+from rich.logging import RichHandler
+
+FORMAT = "%(message)s"
+logging.basicConfig(
+    format=FORMAT, datefmt="[%X]", handlers=[RichHandler(rich_tracebacks=True)]
+)
 
 import numpy as _np
+
+from . import (
+    _skgraph,
+    algorithms,
+    contrib,
+    core,
+    integrators,
+    io,
+    models,
+    network,
+    phase,
+    solvers,
+    topotools,
+    utils,
+    visualization,
+)
+from .utils import Project, Workspace
+
 _np.seterr(divide='ignore', invalid='ignore')
 
 __version__ = utils._get_version()
-
-utils._setup_logger_rich()

openpnm/__version__.py

@@ -1 +1 @@
-__version__ = '3.2.1.dev0'
+__version__ = '3.2.1.dev4'

openpnm/_skgraph/generators/_delaunay.py

@@ -5,8 +5,15 @@
 from openpnm._skgraph.operations import trim_nodes
 
 
-def delaunay(points, shape=[1, 1, 1], reflect=False, trim=True,
-             node_prefix='node', edge_prefix='edge'):
+def delaunay(
+    points,
+    shape=[1, 1, 1],
+    reflect=False,
+    f=1,
+    trim=True,
+    node_prefix='node',
+    edge_prefix='edge',
+):
     r"""
     Generate a network based on Delaunay triangulation of random points
 
@@ -22,6 +29,12 @@ def delaunay(points, shape=[1, 1, 1], reflect=False, trim=True,
         prior to performing the tessellation. These reflected points are
         automatically trimmed.  Enabling this behavior prevents long-range
         connections between surface pores.
+    f : float
+        The fraction of points which should be reflected.  The default is 1 which
+        reflects all the points in the domain, but this can lead to a lot of
+        unnecessary points, so setting to 0.1 or 0.2 helps speed, but risks that
+        the tessellation may not have smooth faces if not enough points are
+        reflected.
     trim : boolean, optional (default = ``True``)
         If ``True`` then any points laying outside the domain are removed. This is
         mostly only useful if ``reflect=True``.
@@ -33,7 +46,7 @@ def delaunay(points, shape=[1, 1, 1], reflect=False, trim=True,
     tri : Delaunay tessellation object
         The Delaunay tessellation object produced by ``scipy.spatial.Delaunay``
     """
-    points = tools.parse_points(points=points, shape=shape, reflect=reflect)
+    points = tools.parse_points(points=points, shape=shape, reflect=reflect, f=f)
     mask = ~np.all(points == 0, axis=0)
     tri = sptl.Delaunay(points=points[:, mask])
     coo = tri_to_am(tri)

openpnm/_skgraph/generators/_voronoi.py

@@ -5,8 +5,16 @@
 from openpnm._skgraph.operations import trim_nodes
 
 
-def voronoi(points, shape=[1, 1, 1], trim=True, reflect=False, relaxation=0,
-            node_prefix='node', edge_prefix='edge'):
+def voronoi(
+    points,
+    shape=[1, 1, 1],
+    trim=True,
+    reflect=False,
+    f=1,
+    relaxation=0,
+    node_prefix='node',
+    edge_prefix='edge',
+):
     r"""
     Generate a network based on a Voronoi tessellation of base points
 
@@ -24,6 +32,12 @@ def voronoi(points, shape=[1, 1, 1], trim=True, reflect=False, relaxation=0,
         If ``True`` then points are reflected across each face of the domain
         prior to performing the tessellation. Enabling this behavior creates
         flat faces on all sizes of the domain.
+    f : float
+        The fraction of points which should be reflected.  The default is 1 which
+        reflects all the points in the domain, but this can lead to a lot of
+        unnecessary points, so setting to 0.1 or 0.2 helps speed, but risks that
+        the tessellation may not have smooth faces if not enough points are
+        reflected.
     relaxation : int, optional (default = 0)
         The number of iterations to use for relaxing the base points. This is
         sometimes called `Lloyd's algorithm
@@ -42,7 +56,7 @@ def voronoi(points, shape=[1, 1, 1], trim=True, reflect=False, relaxation=0,
         The Voronoi tessellation object produced by ``scipy.spatial.Voronoi``
 
     """
-    points = tools.parse_points(points=points, shape=shape, reflect=reflect)
+    points = tools.parse_points(points=points, shape=shape, reflect=reflect, f=f)
     mask = ~np.all(points == 0, axis=0)
     # Perform tessellation
     vor = sptl.Voronoi(points=points[:, mask])

openpnm/_skgraph/generators/_voronoi_delaunay_dual.py

@@ -7,8 +7,17 @@
 from openpnm._skgraph.queries import find_neighbor_nodes
 
 
-def voronoi_delaunay_dual(points, shape, trim=True, reflect=True, relaxation=0,
-                          node_prefix='node', edge_prefix='edge'):
+def voronoi_delaunay_dual(
+    points,
+    shape,
+    trim=True,
+    reflect=True,
+    f=1,
+    relaxation=0,
+    node_prefix='node',
+    edge_prefix='edge',
+    return_tri=False,
+):
     r"""
     Generate a dual Voronoi-Delaunay network from given base points
 
@@ -25,6 +34,12 @@
     reflect : bool, optionl
         If ``True`` (Default) then points are reflected across each face of the
         domain prior to performing the tessellation.
+    f : float
+        The fraction of points which should be reflected.  The default is 1 which
+        reflects all the points in the domain, but this can lead to a lot of
+        unnecessary points, so setting to 0.1 or 0.2 helps speed, but risks that
+        the tessellation may not have smooth faces if not enough points are
+        reflected.
     relaxation : int, optional (default = 0)
         The number of iterations to use for relaxing the base points. This is
         sometimes called `Lloyd's algorithm
@@ -42,11 +57,17 @@
     vor : Voronoi object
         The Voronoi tessellation object produced by ``scipy.spatial.Voronoi``
     tri : Delaunay object
-        The Delaunay triangulation object produced ``scipy.spatial.Delaunay``
+        The Delaunay triangulation object produced by ``scipy.spatial.Delaunay``
 
     """
     # Generate a set of base points if scalar was given
-    points = tools.parse_points(points=points, shape=shape, reflect=reflect)
+    points = tools.parse_points(
+        points=points,
+        shape=shape,
+        reflect=reflect,
+        f=1,
+    )
+
     # Generate mask to remove any dims with all 0's
     mask = ~np.all(points == 0, axis=0)
 
@@ -55,45 +76,36 @@
     for _ in range(relaxation):
         points = tools.lloyd_relaxation(vor, mode='fast')
         vor = sptl.Voronoi(points=points[:, mask])
-    tri = sptl.Delaunay(points=points[:, mask])
-
-    # Combine points
-    pts_all = np.vstack((vor.points, vor.vertices))
-    Nall = np.shape(pts_all)[0]
-
-    # Create adjacency matrix in lil format for quick construction
-    am = sprs.lil_matrix((Nall, Nall))
-    for ridge in vor.ridge_dict.keys():
-        # Make Delaunay-to-Delaunay connections
-        for i in ridge:
-            am.rows[i].extend([ridge[0], ridge[1]])
-        # Get Voronoi vertices for current ridge
-        row = vor.ridge_dict[ridge].copy()
-        # Index Voronoi vertex numbers by number of Delaunay points
-        row = [i + vor.npoints for i in row if i > -1]
-        # Make Voronoi-to-Delaunay connections
-        for i in ridge:
-            am.rows[i].extend(row)
-        # Make Voronoi-to-Voronoi connections
-        row.append(row[0])
-        for i in range(len(row)-1):
-            am.rows[row[i]].append(row[i+1])
-
-    # Finalize adjacency matrix by assigning data values
-    am.data = am.rows  # Values don't matter, only shape, so use 'rows'
-    # Convert to COO format for direct acces to row and col
-    am = am.tocoo()
-    # Extract rows and cols
-    conns = np.vstack((am.row, am.col)).T
 
-    # Convert to sanitized adjacency matrix
+    # Collect delaunay edges
+    conns_del = vor.ridge_points
+    # Deal with voronoi edges
+    v = vor.ridge_vertices.copy()  # Assuming 'ridge' means facet between regions
+    # Add row [0] to close the facet on itself, add -1 to break connection to
+    # next facet in list as connections with -1 get deleted later
+    _ = [row.extend([row[0], -1]) for row in v]
+    v = np.hstack(v)
+    conns_vor = np.vstack((v[:-1], v[1:])).T
+    mask = np.any(conns_vor < 0, axis=1)
+    conns_vor = conns_vor[~mask] + vor.npoints
+    # Finally, get interconnecting edges
+    idx = [vor.regions[vor.point_region[i]] for i in range(0, len(vor.regions)-1)]
+    conns_inter = [([i]*len(idx[i]), idx[i]) for i in range(0, len(idx))]
+    conns_inter = np.hstack(conns_inter).astype(int).T
+    mask = np.any(conns_inter < 0, axis=1)
+    conns_inter = conns_inter[~mask, :] + np.array([0, vor.npoints], dtype=int)
+    conns = np.vstack((conns_del, conns_vor, conns_inter))
+
+    # Tidy up
     am = conns_to_am(conns)
-    # Finally, retreive conns back from am
     conns = np.vstack((am.row, am.col)).T
 
-    # Convert coords to 3D if necessary
+    # Combine delaunay and voronoi points
+    pts_all = np.vstack((vor.points, vor.vertices))
     # Rounding is crucial since some voronoi verts endup outside domain
     pts_all = np.around(pts_all, decimals=10)
+    # Convert coords to 3D if necessary
+    mask = ~np.all(points == 0, axis=0)
     if mask.sum() < 3:
         coords = np.zeros([pts_all.shape[0], 3], dtype=float)
         coords[:, mask] = pts_all
@@ -142,4 +154,35 @@
             inds = np.where(trim)[0]
             network = trim_nodes(network=network, inds=inds)
 
+    if return_tri:
+        tri = sptl.Delaunay(points=points[:, mask])
+    else:
+        tri = None
     return network, vor, tri
+
+
+if __name__ == "__main__":
+    import openpnm as op
+    pn, vor, tri = voronoi_delaunay_dual(
+        points=500,
+        shape=[1, 1, 0],
+        trim=True,
+        reflect=True,
+        f=0.2,
+        relaxation=0,
+        node_prefix='pore',
+        edge_prefix='throat',
+    )
+    net = op.network.Network()
+    net.update(pn)
+    h = None
+    h = op.visualization.plot_connections(
+        net, throats=pn['throat.delaunay'], c='b', ax=h)
+    h = op.visualization.plot_connections(
+        net, throats=pn['throat.voronoi'], c='g', ax=h)
+    h = op.visualization.plot_connections(
+        net, throats=pn['throat.interconnect'], c='r', ax=h)
+    h = op.visualization.plot_coordinates(
+        net, pores=pn['pore.voronoi'], c='g', markersize=150, ax=h)
+    h = op.visualization.plot_coordinates(
+        net, pores=pn['pore.delaunay'], c='b', markersize=150, ax=h)