Stitching parts of same network

[aaornsc]

Hi, I've been using OpenPNM for a few months by now. The new tutorials and examples are really helpful.

I have some CT image data, I extracted the network and imported it to OpenPNM. The network is connected at first but then i trimmed the network to a new domain and the network is not connected anymore.

This is the trimmed network:

I want the following areas to be connected: 1 to 2 and 2 to 3.

In the manipulation examples there are two networks as a starting point when stitching, so i know i have to divide or make subsets of my trimmed network to have subnetworks and then stitch them. My question is how do i label each pore and throat? Any function in topotools is helpful in my case? Or do i have to go pore by pore and throat by throat and label them? Also i don't want any pore or throat outside the areas 1,2 and 3 to be part of the new connected network.

Thanks for the help!!

[jgostick]

So, your act of trimming the network is breaking the connectivity. It's possible to restitch these together but this is not likely to be representative of your original material.

Assuming you know the above point, it "might" be possible to stitch these together but OpenPNM does not offer a built in function for this. We have something similar, topotools.stitch, but this works on separate networks, while your case they are already merged into one. I think that we should add something like this to openpnm. I will attempt to write a function for this soon. It shouldn't be too bad, we have the logic in the stitch method already.

[jgostick]

OK, I have written up something that does the job. I think it will be a useful addition to openpnm, so didn't mind investing a bit of time on it. Below is a script that illlustrates the idea, and I'll probably add the function to topotools soon, but for now you can try this. You have to identify the pores in each of the isolated clusters yourself though:

import openpnm as op
import numpy as np
from openpnm.topotools import extend

def stitch_pores(network, pores1, pores2, mode='gabriel'):
    from openpnm.network import Delaunay, Gabriel
    pores1 = network._parse_indices(pores1)
    pores2 = network._parse_indices(pores2)
    C1 = network.coords[pores1, :]
    C2 = network.coords[pores2, :]
    crds = np.vstack((C1, C2))
    if mode == 'delaunay':
        net = Delaunay(points=crds, settings={'trim': False})
    if mode == 'gabriel':
        net = Gabriel(points=crds, settings={'trim': False})
    net.set_label(pores=range(len(pores1)), label='pore.one')
    net.set_label(pores=range(len(pores2)), label='pore.two')
    Ts = net.find_neighbor_throats(pores=net.pores('one'), mode='xor')
    conns = net.conns[Ts]
    mapped_conns = np.vstack((pores1[conns[:, 0]],
                              pores2[conns[:, 1] - len(pores1)])).T
    mapped_conns = np.sort(mapped_conns, axis=1)
    extend(network=network, conns=mapped_conns, labels='stitched')

pn = op.network.Delaunay(num_points=500)
Ps = (pn.coords[:, 0] > 0.3) * (pn.coords[:, 0] < 0.6)
op.topotools.trim(network=pn, pores=Ps)
pores1 = pn.coords[:, 0] < 0.3
pores2 = pn.coords[:, 0] > 0.6
h = pn.check_network_health()
stitch_pores(network=pn, pores1=pores1, pores2=pores2, mode='gabriel')
fig = op.topotools.plot_coordinates(pn, markersize=50, c='g')
fig = op.topotools.plot_connections(pn, throats=pn.throats('stitched', mode='not'), fig=fig)
fig = op.topotools.plot_connections(pn, throats=pn.throats('stitched'), c='r', fig=fig)

[jgostick]

[aaornsc]

Wow, I really appreciate the help, I'll let you know as soon as I use the function on my network. I still have to figure out how to identify the pores i want to connect, but since they are less than 250 i can just work with the visualization on Paraview and start from there. Thanks for the comments above too.