misc.grid.unstruct

  1"""\
  2Convert cornerpoint grids to unstructured grids
  3
  4:example:
  5    import pyresito.grid.unstruct as us
  6    import pyresito.io.grdecl as grdecl
  7
  8    g = grdecl.read('~/proj/cmgtools/bld/overlap.grdecl')
  9"""
 10# pylint: disable=too-few-public-methods, multiple-statements
 11import numpy as np
 12
 13
 14class Ridge (object):
 15    """A ridge consists of two points, anchored in each their pillar. We only
 16    need to store the z-values, because the x- and y- values are determined by
 17    the pillar themselves.
 18    """
 19    __slots__ = ['left', 'right']
 20
 21    def __init__(self, left, right):
 22        self.left = left
 23        self.right = right
 24
 25    def is_not_below(self, other):
 26        """Weak ordering of ridges based on vertical placement.
 27
 28        :param other:  Ridge to be compared to this object.
 29        :type other:   :class:`Ridge`
 30        :returns:      True if no point on self is below any on the other,
 31                       None if the ridges cross, and False if there is a point
 32                       on the other ridge that is above any on self.
 33        :rtype:        Boolean
 34        """
 35        # test each side separately. if self is at the same level as other,
 36        # this should count positively towards the test (i.e. it is regarded
 37        # as "at least as high"), so use less-or-equal.
 38        left_above = other.left <= self.left
 39        right_above = other.right <= self.right
 40
 41        # if both sides are at least as high, then self is regarded as at
 42        # least as high as other, and vice versa. if one side is lower and
 43        # one side is higher, then the ridges cross.
 44        if left_above:
 45            return True if right_above else None
 46        else:
 47            return None if right_above else False
 48
 49
 50class Face (object):
 51    """A (vertical) face consists of two ridges, because all of the faces in a
 52    hexahedron can be seen as (possibly degenerate) quadrilaterals.
 53    """
 54    __slots__ = ['top', 'btm']
 55
 56    def __init__(self, top, btm):
 57        self.top = top
 58        self.btm = btm
 59
 60    def is_above(self, other):
 61        """Weak ordering of faces based on vertical placement.
 62
 63        :param other:   Face to be compared to this object.
 64        :type other:    :class:`Face`
 65        :returns:       True if all points in face self is above all points in
 66                        face other, False otherwise
 67        :rtype:         Boolean
 68        """
 69        # if the bottom of self is aligned with the top of other, then the
 70        # face itself is considered to be above. since the ridge test also
 71        # has an indeterminate result, we test explicitly like this
 72        return True if self.btm.is_not_below(other.top) else False
 73
 74
 75def conv(grid):
 76    """Convert a cornerpoint grid to an unstructured grid.
 77
 78    :param grid:  Cornerpoint grid to be converted.
 79    :type grid:   dict, with 'COORD', 'ZCORN', 'ACTNUM'
 80    :returns:     Unstructured grid
 81    :rtype:       dict
 82    """
 83    # extract the properties of interest from the cornerpoint grid
 84    zcorn = grid['ZCORN']
 85    actnum = grid['ACTNUM']
 86    ni, nj, nk = grid['DIMENS']
 87
 88    # zcorn has now dimensionality (k, b, j, f, i, r) and actnum is (k, j, i);
 89    # permute the cubes to get the (b, k)-dimensions varying quickest, to avoid
 90    # page faults when we move along pillars/columns
 91    zcorn = np.transpose(zcorn,  axes=[2, 3, 4, 5, 1, 0])
 92    actnum = np.transpose(actnum, axes=[1, 2, 0])
 93
 94    # memory allocation: number of unique cornerpoints along each pillar, and
 95    # the index of each cornerpoint into the global list of vertices
 96    num_cp = np.empty((nj + 1, ni + 1), dtype=np.int32)
 97    ndx_cp = np.empty((nk, 2, nj, 2, ni, 2), dtype=np.int32)
 98
 99    # each pillar is connected to at most 2*2 columns (front/back, right/left),
100    # and each column has at most 2*nk (one top and one bottom) corners
101    corn_z = np.empty((2, 2, 2, nk), dtype=np.float32)
102    corn_i = np.empty((2, 2, 2, nk), dtype=np.int32)
103    corn_j = np.empty((2, 2, 2, nk), dtype=np.int32)
104    corn_a = np.empty((2, 2, 2, nk), dtype=np.bool)
105
106    # get all unique points that are hinged to a certain pillar (p, q)
107    for q, p in np.ndindex((nj + 1, ni + 1)):
108        # reset number of corners found for this column
109        num_corn_z = 0
110
111        for f, r in np.ndindex((2, 2)):  # front/back, right/left
112            # calculate the cell index of the column at this position
113            j = q - f
114            i = p - r
115
116            for b in range(2):  # bottom/top
117
118                # copy depth values for this corner; notice that we have
119                # pivoted the zcorn matrix so that values going upwards are in
120                # last dim.
121                corn_z[f, r, b, :] = zcorn[j, f, i, r, b, :]
122
123                # same for active numbers, but this is reused for top/bottom;
124                # if the cell is inactive, then both top and bottom point are.
125                corn_a[f, r, b, :] = actnum[j, i, :]
126
127                # save the original indices into these auxiliary arrays so that
128                # we can figure out where each point came from after they are
129                # sorted
130                corn_i[f, r, b] = i
131                corn_j[f, r, b] = j