<commit_before>import scipy as sp
import numpy as np
import scipy.spatial as sptl
import matplotlib.pyplot as plt
from transforms3d import _gohlketransforms as tr
from scipy import ndimage
import math
from skimage.morphology import convex_hull_image
from skimage.measure import regionprops
from openpnm import topotools
from openpnm.network import DelaunayVoronoiDual
from openpnm.utils import logging, Project
import openpnm.models.geometry as gm
from openpnm.geometry import GenericGeometry
from openpnm.utils.misc import unique_list
import openpnm.utils.vertexops as vo
from scipy.stats import itemfreq
logger = logging.getLogger(__name__)


class VoronoiFibers(Project):
    r"""
    Resembles a fibrous paper or mat with straight intersecting fibers.

    Two geometries are created: DelaunayGeometry defines the pore space
    with pores connected by a Delaunay tesselation and VoronoiGeometry defines
    the fiber space with fibers forming the edges of the Voronoi diagram.
    The two geometries are complimentary and can be accessed individually
    via the project associated with this network object.

    This class is a subclass of the DelaunayVoronoiDual network, and many of
    the parameters behave in the exact same way. However, currently a
    rectangular shape is the only one supported. Additionally, an image is
    created for calculating size data and two parameters relate to this:
    fiber_rad and resolution as detailed below.

    Parameters
    ----------
    num_points : integer
        The number of random base points to distribute inside the domain.
        These points will become connected by the Delaunay triangulation.  The
        points will be generated by calling ``generate_base_points`` in
        ``topotools``.

    points : array_like (num_points x 3)
        A list of coordinates for pre-generated points, typically produced
        using ``generate_base_points`` in topotools.  Note that base points
        should extend beyond the domain so that degenerate Voronoi points
        can be trimmed.  Note: the spherical and cylindrical options
        cannot be used here.

    shape : array_like
        The size and shape of the domain using for generating and trimming
        excess points. It's treated as the outer corner of rectangle [x, y, z]
        whose opposite corner lies at [0, 0, 0].

        By default, a domain size of [1, 1, 1] is used.

    fiber_rad: float
        fiber radius to apply to Voronoi edges when calculating pore and throat
        sizes

    resolution : boolean
        Determines the size of each voxel in the image. Care should be made to
        appropriately set the resolution based on the fiber_radius and the
        shape of the domain so as to remain within memory constraints.

    References
    ----------
    This approach to modeling fibrous materials was first presented by
    Thompson [1] for simulating fluid imbibition in sorbent paper products.
    Gostick [2], and Tranter et al.[3, 4] have subsequently used it to model
    electrodes in fuel cells.

    [1] <NAME>, AlChE J., 48, 1369 (2002)
    [2] <NAME>, Journal of the Electrochemical Society 2013, 160, F731.
    [3] <NAME>. Tranter et al. Fuel Cells, 2016, 16, 4, 504-515
    [4] <NAME>. Tranter et al. Transport in Porous Media, 2018, 121, 3, 597-620

    Examples
    --------
    Points will be automatically generated if none are given:

    >>> import openpnm as op
    >>> ws = op.Workspace()
    >>> ws.clear()
    >>> prj = op.materials.VoronoiFibers(num_points=50,
    ...                                  shape=[1e-4, 1e-4, 1e-4],
    ...                                  fiber_rad=5e-6,
    ...                                  resolution=1e-6)
    """

    def __init__(self, num_points=None, points=None, shape=[1, 1, 1],
                 fiber_rad=None, resolution=1e-2, name=None, **kwargs):
        super().__init__(name=name)
        shape = np.array(shape)
        if (len(shape) != 3) or np.any(shape == 0):
            raise Exception('Only 3D, rectangular shapes are supported')
        if fiber_rad is None:
            logger.exception(msg='Please initialize class with a fiber_rad')

        net = DelaunayVoronoiDual(project=self,
                                  num_points=num_points,
                                  points=points,
                                  shape=shape,
                                  name=self.name+'_net',
                                  **kwargs)
        net.fiber_rad = fiber_rad
        net.resolution = resolution
        DelaunayGeometry(project=self,
                         network=net,
                         pores=net.pores('delaunay'),
                         throats=net.throats('delaunay'),
                         name=self.name+'_del')
        VoronoiGeometry(project=self,
                        network=net,
                        pores=net.pores('voronoi'),
                        throats=net.throats('voronoi'),
                        name=self.name+'_vor')


class DelaunayGeometry(GenericGeometry):
    r"""
    Subclass of GenericGeometry for the pores connected by throats formed from
    the Delaunay tesselation.

    Parameters
    ----------
    name : string
        A unique name for the network
    """

    def __init__(self, network, **kwargs):
        super().__init__(network=network, **kwargs)
        # Set all the required models
        vertices = network.find_pore_hulls()
        p_coords = np.array([network['pore.coords'][list(p)] for p in vertices],
                            dtype=object)
        self['pore.vertices'] = p_coords
        vertices = network.find_throat_facets()
        t_coords = np.array([network['pore.coords'][list(t)] for t in vertices],
                            dtype=object)
        self['throat.vertices'] = t_coords

        self.in_hull_volume()
        self['throat.normal'] = self._t_normals()
        self['throat.centroid'] = self._centroids(verts=t_coords)
        self['pore.centroid'] = self._centroids(verts=p_coords)
        (self['pore.indiameter'],
         self['pore.incenter']) = self._indiameter_from_fibers()
        self._throat_props()
        topotools.trim_occluded_throats(network=network, mask=self.name)
        self['throat.volume'] = sp.zeros(1, dtype=float)
        self['throat.length'] = sp.ones(1, dtype=float)
        self['throat.length'] *= self.network.fiber_rad*2
        cen_lens = self._throat_c2c()
        self['throat.c2c'] = np.sum(cen_lens, axis=1)
        cen_lens[cen_lens <= 0.0] = 1e-12
        self['throat.centroid_lengths.pore1'] = cen_lens[:, 0]
        self['throat.centroid_lengths.throat'] = cen_lens[:, 1]
        self['throat.centroid_lengths.pore2'] = cen_lens[:, 2]
        # Configurable Models
        self.add_model(propname='throat.equivalent_area',
                       model=gm.throat_equivalent_area.boundary)
        self.add_model(propname='throat.shape_factor',
                       model=gm.throat_shape_factor.compactness)
        self.add_model(propname='pore.diameter',
                       model=gm.pore_size.equivalent_diameter)
        self.add_model(propname='pore.area',
                       model=gm.pore_area.sphere,
                       pore_diameter='pore.diameter')
        self.add_model(propname='throat.surface_area',
                       model=gm.throat_surface_area.extrusion)
        self.regenerate_models()

    def _t_normals(self):
        r"""
        Update the throat normals from the voronoi vertices
        """
        verts = self['throat.vertices']
        value = sp.zeros([len(verts), 3])
        for i in range(len(verts)):
            if len(sp.unique(verts[i][:, 0])) == 1:
                verts_2d = sp.vstack((verts[i][:, 1], verts[i][:, 2])).T
            elif len(sp.unique(verts[i][:, 1])) == 1:
                verts_2d = sp.vstack((verts[i][:, 0], verts[i][:, 2])).T
            else:
                verts_2d = sp.vstack((verts[i][:, 0], verts[i][:, 1])).T
            hull = sptl.ConvexHull(verts_2d, qhull_options='QJ Pp')
            sorted_verts = verts[i][hull.vertices].astype(float)
            v1 = sorted_verts[-1]-sorted_verts[0]
            v2 = sorted_verts[1]-sorted_verts[0]
            value[i] = tr.unit_vector(sp.cross(v1, v2))
        return value

    def _centroids(self, verts):
        r'''
        Function to calculate the centroid as the mean of a set of vertices.
        Used for pore and throat.
        '''
        value = sp.zeros([len(verts), 3])
        for i, i_verts in enumerate(verts):
            value[i] = np.mean(i_verts, axis=0)<commit_msg>Fix typo and rearrange models<commit_after>import scipy as sp
import numpy as np
import scipy.spatial as sptl
import matplotlib.pyplot as plt
from transforms3d import _gohlketransforms as tr
from scipy import ndimage
import math
from skimage.morphology import convex_hull_image
from skimage.measure import regionprops
from openpnm import topotools
from openpnm.network import DelaunayVoronoiDual
from openpnm.utils import logging, Project
import openpnm.models.geometry as gm
from openpnm.geometry import GenericGeometry
from openpnm.utils.misc import unique_list
import openpnm.utils.vertexops as vo
from scipy.stats import itemfreq
logger = logging.getLogger(__name__)


class VoronoiFibers(Project):
    r"""
    Resembles a fibrous paper or mat with straight intersecting fibers.

    Two geometries are created: DelaunayGeometry defines the pore space
    with pores connected by a Delaunay tesselation and VoronoiGeometry defines
    the fiber space with fibers forming the edges of the Voronoi diagram.
    The two geometries are complimentary and can be accessed individually
    via the project associated with this network object.

    This class is a subclass of the DelaunayVoronoiDual network, and many of
    the parameters behave in the exact same way. However, currently a
    rectangular shape is the only one supported. Additionally, an image is
    created for calculating size data and two parameters relate to this:
    fiber_rad and resolution as detailed below.

    Parameters
    ----------
    num_points : integer
        The number of random base points to distribute inside the domain.
        These points will become connected by the Delaunay triangulation.  The
        points will be generated by calling ``generate_base_points`` in
        ``topotools``.

    points : array_like (num_points x 3)
        A list of coordinates for pre-generated points, typically produced
        using ``generate_base_points`` in topotools.  Note that base points
        should extend beyond the domain so that degenerate Voronoi points
        can be trimmed.  Note: the spherical and cylindrical options
        cannot be used here.

    shape : array_like
        The size and shape of the domain using for generating and trimming
        excess points. It's treated as the outer corner of rectangle [x, y, z]
        whose opposite corner lies at [0, 0, 0].

        By default, a domain size of [1, 1, 1] is used.

    fiber_rad: float
        fiber radius to apply to Voronoi edges when calculating pore and throat
        sizes

    resolution : boolean
        Determines the size of each voxel in the image. Care should be made to
        appropriately set the resolution based on the fiber_radius and the
        shape of the domain so as to remain within memory constraints.

    References
    ----------
    This approach to modeling fibrous materials was first presented by
    Thompson [1] for simulating fluid imbibition in sorbent paper products.
    Gostick [2], and Tranter et al.[3, 4] have subsequently used it to model
    electrodes in fuel cells.

    [1] <NAME>, AlChE J., 48, 1369 (2002)
    [2] <NAME>, Journal of the Electrochemical Society 2013, 160, F731.
    [3] <NAME>. Tranter et al. Fuel Cells, 2016, 16, 4, 504-515
    [4] <NAME> et al. Transport in Porous Media, 2018, 121, 3, 597-620

    Examples
    --------
    Points will be automatically generated if none are given:

    >>> import openpnm as op
    >>> ws = op.Workspace()
    >>> ws.clear()
    >>> prj = op.materials.VoronoiFibers(num_points=50,
    ...                                  shape=[1e-4, 1e-4, 1e-4],
    ...                                  fiber_rad=5e-6,
    ...                                  resolution=1e-6)
    """

    def __init__(self, num_points=None, points=None, shape=[1, 1, 1],
                 fiber_rad=None, resolution=1e-2, name=None, **kwargs):
        super().__init__(name=name)
        shape = np.array(shape)
        if (len(shape) != 3) or np.any(shape == 0):
            raise Exception('Only 3D, rectangular shapes are supported')
        if fiber_rad is None:
            logger.exception(msg='Please initialize class with a fiber_rad')

        net = DelaunayVoronoiDual(project=self,
                                  num_points=num_points,
                                  points=points,
                                  shape=shape,
                                  name=self.name+'_net',
                                  **kwargs)
        net.fiber_rad = fiber_rad
        net.resolution = resolution
        DelaunayGeometry(project=self,
                         network=net,
                         pores=net.pores('delaunay'),
                         throats=net.throats('delaunay'),
                         name=self.name+'_del')
        VoronoiGeometry(project=self,
                        network=net,
                        pores=net.pores('voronoi'),
                        throats=net.throats('voronoi'),
                        name=self.name+'_vor')


class DelaunayGeometry(GenericGeometry):
    r"""
    Subclass of GenericGeometry for the pores connected by throats formed from
    the Delaunay tesselation.

    Parameters
    ----------
    name : string
        A unique name for the network
    """

    def __init__(self, network, **kwargs):
        super().__init__(network=network, **kwargs)
        # Set all the required models
        vertices = network.find_pore_hulls()
        p_coords = np.array([network['pore.coords'][list(p)] for p in vertices],
                            dtype=object)
        self['pore.vertices'] = p_coords
        vertices = network.find_throat_facets()
        t_coords = np.array([network['pore.coords'][list(t)] for t in vertices],
                            dtype=object)
        self['throat.vertices'] = t_coords

        self.in_hull_volume()
        self['throat.normal'] = self._t_normals()
        self['throat.centroid'] = self._centroids(verts=t_coords)
        self['pore.centroid'] = self._centroids(verts=p_coords)
        (self['pore.indiameter'],
         self['pore.incenter']) = self._indiameter_from_fibers()
        self._throat_props()
        topotools.trim_occluded_throats(network=network, mask=self.name)
        self['throat.volume'] = sp.zeros(1, dtype=float)
        self['throat.length'] = sp.ones(1, dtype=float)
        self['throat.length'] *= self.network.fiber_rad*2
        cen_lens = self._throat_c2c()
        self['throat.c2c'] = np.sum(cen_lens, axis=1)
        cen_lens[cen_lens <= 0.0] = 1e-12
        self['throat.conduit_lengths.pore1'] = cen_lens[:, 0]
        self['throat.conduit_lengths.throat'] = cen_lens[:, 1]
        self['throat.conduit_lengths.pore2'] = cen_lens[:, 2]
        # Configurable Models
        self.add_model(propname='throat.shape_factor',
                       model=gm.throat_shape_factor.compactness)
        self.add_model(propname='pore.diameter',
                       model=gm.pore_size.equivalent_diameter)
        self.add_model(propname='pore.area',
                       model=gm.pore_area.sphere,
                       pore_diameter='pore.diameter')
        self.add_model(propname='throat.surface_area',
                       model=gm.throat_surface_area.extrusion)
        self.add_model(propname='throat.equivalent_area',
                       model=gm.throat_equivalent_area.boundary)
        self.regenerate_models()

    def _t_normals(self):
        r"""
        Update the throat normals from the voronoi vertices
        """
        verts = self['throat.vertices']
        value = sp.zeros([len(verts), 3])
        for i in range(len(verts)):
            if len(sp.unique(verts[i][:, 0])) == 1:
                verts_2d = sp.vstack((verts[i][:, 1], verts[i][:, 2])).T
            elif len(sp.unique(verts[i][:, 1])) == 1:
                verts_2d = sp.vstack((verts[i][:, 0], verts[i][:, 2])).T
            else:
                verts_2d = sp.vstack((verts[i][:, 0], verts[i][:, 1])).T
            hull = sptl.ConvexHull(verts_2d, qhull_options='QJ Pp')
            sorted_verts = verts[i][hull.vertices].astype(float)
            v1 = sorted_verts[-1]-sorted_verts[0]
            v2 = sorted_verts[1]-sorted_verts[0]
            value[i] = tr.unit_vector(sp.cross(v1, v2))
        return value

    def _centroids(self, verts):
        r'''
        Function to calculate the centroid as the mean of a set of vertices.
        Used for pore and throat.
        '''
        value = sp.zeros([len(verts), 3])
        for i, i_verts in enumerate(verts):
            value[i] = np.mean(i_verts, axis=0)
{'fraction_non_alphanumeric': 0.09750580046403712, 'fraction_numerical': 0.014385150812064965, 'mean_word_length': 2.2413987591652567, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2117908', 'n_tokens_mistral': 5389, 'n_tokens_neox': 4981, 'n_words': 1452}