import warnings
import scipy as sp
import numpy as np
import porespy as ps
import openpnm as op
import matplotlib.pyplot as plt
ws = op.Workspace()
ws.settings["loglevel"] = 40
warnings.filterwarnings('ignore')
%matplotlib inline
np.random.seed(10)

# NBVAL_IGNORE_OUTPUT
im = ps.generators.overlapping_spheres(shape=[200, 200, 200], r=10, porosity=0.5, maxiter=0)
plt.imshow(im[:, :, 50]);

"""
Explanation: Working with Extracted Networks
Extracting a pore network using PoreSpy and loading into OpenPNM
End of explanation
"""


eps = ps.metrics.porosity(im)
print(f"Porosity: {eps*100:.1f}%")

"""
Explanation: Let's check out the porosity of the generated image!
End of explanation
"""


# NBVAL_IGNORE_OUTPUT
im_3d = ps.visualization.show_3D(im)
plt.imshow(im_3d, cmap=plt.cm.magma);

# NBVAL_IGNORE_OUTPUT
snow = ps.networks.snow2(im, boundary_width=[[0, 3], 0, 0])

"""
Explanation: Let's visualize the image using porespy's 3D visualizer: (this might take several seconds)
End of explanation
"""


# NBVAL_IGNORE_OUTPUT
proj = op.io.PoreSpy.import_data(snow.network)
print(proj)

"""
Explanation: OpenPNM has an IO class specifically for importing the output from PoreSpy.  The import_data method can either accept a handle to a dictionary (as output from the snow algorithm above), or it can accept a filename to a saved dctionary (saved using Python's pickle library).  All IO methods in OpenPNM return a project which is a list, in this case containing a network and a geometry object.
End of explanation
"""


# NBVAL_IGNORE_OUTPUT
net = proj[0]
geo = proj[1]
print(net)

"""
Explanation: We can unpack the network and geometry objects from the project using the indices in the list as follows:
End of explanation
"""


# NBVAL_IGNORE_OUTPUT
print(geo)

"""
Explanation: It is important to note that the net object only has topological information and labels.  The geo object was created by the openpnm.io.PoreSpy import class to extract all geometric information from the supplied snow dict and put in on a geometry object.  We can print geo to confirm:
End of explanation
"""


# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(1, 1, figsize=[8, 8])
ax = op.topotools.plot_connections(network=net, alpha=0.8, color='grey', ax=ax)
ax = op.topotools.plot_coordinates(network=net, ax=ax, color='b', markersize=50)
ax = op.topotools.plot_coordinates(network=net, pores=net.pores('xmax'), ax=ax, color='r', markersize=50)

"""
Explanation: Now let's plot things to see what we have:
End of explanation
"""


m = np.meshgrid(range(50, 195, 10), range(50, 195, 10))
m = np.vstack([-10*np.ones_like(m[0].flatten()), m[0].flatten(), m[1].flatten()]).T

"""
Explanation: This looks pretty good, but it only has boundary pores on the right face, indicated by the red dots.  When we ran the snow algorithm we specifically told it to only put boundary pores the "right".  We could have added them to all faces during the extraction, but for the sake of demonstration we can add them after the fact, although the result is slightly different, as you'll see.  
We'll use the find_surface_pores function in the topotools module.  This function applies a Delaunay tessellation between the interior pores and some fictitious "marker" nodes.  Only pores that are on the surface will be connected to these marker nodes.  To get the best result from the find_surface_pores function, it's a good idea to supply your own markers, so let's make a 2D plane of points, positioned outside the left face of the domain:
End of explanation
"""


op.topotools.find_surface_pores(network=net, markers=m, label='left')

"""
Explanation: Now we pass these points in as markers to the find_surface_pores function:
End of explanation
"""


op.topotools.clone_pores(network=net, pores=net.pores('left'), labels='left_boundary')
net['pore.coords'][net.pores('left_boundary')] *= [0, 1, 1]

"""
Explanation: Lastly we want to "clone" these pores and translate them to domain edge:
End of explanation
"""


# NBVAL_IGNORE_OUTPUT
Ps = net.pores('xmax')
net.set_label('right_boundary', pores=Ps)
Ts = net.find_neighbor_throats(pores=net.pores('right_boundary'), mode='or')
net.set_label('right_boundary', throats=Ts)

fig, ax = plt.subplots(1, 1, figsize=[5, 5])
ax = op.topotools.plot_coordinates(network=net, color='g', alpha=0.2, ax=ax)
ax = op.topotools.plot_connections(network=net, throats=net.throats('right_boundary'), color='r', ax=ax)
ax = op.topotools.plot_connections(network=net, throats=net.throats('left_boundary'), color='b', ax=ax)

"""
Explanation: Now let's inspect the result using the quick plotting tools in the topotools module.  First we'll add a new label called 'right_boundary' to match the 'left_boundary' we added above, then we'll plot the throats that connect to ther 'right_boundary' or 'left_boundary':
End of explanation
"""


