new logic for vor-to-am

openpnm/_skgraph/generators/_voronoi_delaunay_dual.py

@@ -55,45 +55,38 @@ def voronoi_delaunay_dual(points, shape, trim=True, reflect=True, relaxation=0,
     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
+    # tri = sptl.Delaunay(points=points[:, mask])
+    tri = None
+
+    # 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
+    _ = [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
@@ -143,3 +136,57 @@ def voronoi_delaunay_dual(points, shape, trim=True, reflect=True, relaxation=0,
             network = trim_nodes(network=network, inds=inds)
 
     return network, vor, tri
+
+
+if __name__ == "__main__":
+    import openpnm as op
+    pn, vor, tri = voronoi_delaunay_dual(
+        points=10000,
+        shape=[1, 1, 1],
+        trim=True,
+        reflect=True,
+        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)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+