"""Classes and routines for handling model grids."""
import warnings
import numpy as np
import vtk
from vtkmodules.util import numpy_support
from .decorators import cached_property, apply_to_each_input, state_check
from .base_spatial import SpatialComponent
from .grid_utils import get_top_z_coords, numba_get_volumes, numba_get_xyz, get_connectivity_matrix
from .utils import rolling_window, mk_vtk_id_list, get_single_path
from .parse_utils import read_ecl_bin
[docs]
class Grid(SpatialComponent):
"""Basic grid class."""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self._vtk_grid = None
self._cell_id_d = None
self._vtk_grid_params = {'use_only_active': True, 'cell_size': None, 'scaling': True}
self._vtk_grid_scales = [1, 1, 1]
self._vtk_locator = None
if 'MAPAXES' not in self:
setattr(self, 'MAPAXES', np.array([0, 1, 0, 0, 1, 0]))
if 'ACTNUM' not in self and 'DIMENS' in self:
self.actnum = np.ones(self.dimens, dtype=bool)
@property
def origin(self):
"""Grid axes origin relative to the map coordinates."""
return np.array([self.mapaxes[2], self.mapaxes[3], self.tops])
@cached_property
def bounding_box(self):
"""Pair of diagonal corner points for grid's bounding box."""
raise NotImplementedError()
@cached_property
def cell_centroids(self):
"""Centroids of cells."""
raise NotImplementedError()
@cached_property
def cell_volumes(self):
"""Volumes of cells."""
grid = specify_grid(self)
return grid.cell_volumes
@cached_property
def xyz(self):
"""Cell coordinates."""
raise NotImplementedError()
@property
def ex(self):
"""Unit vector along grid X axis."""
ex = np.array([self.mapaxes[-2] - self.mapaxes[2], self.mapaxes[-1] - self.mapaxes[3]])
return ex / np.linalg.norm(ex)
@property
def ey(self):
"""Unit vector along grid Y axis."""
ey = np.array([self.mapaxes[0] - self.mapaxes[2], self.mapaxes[1] - self.mapaxes[3]])
return ey / np.linalg.norm(ey)
def _load_ecl_binary(self, path_to_results, attrs, basename, logger=None,
**kwargs):
_ = kwargs
path = get_single_path(path_to_results, basename + '.EGRID', logger)
if path is None:
return
attrs_tmp = attrs + ['GRIDHEAD'] if 'DIMENS' in attrs else attrs
sections = read_ecl_bin(path, attrs_tmp, logger=logger)
if 'DIMENS' in attrs:
setattr(self, 'DIMENS', sections['GRIDHEAD'][1:4])
for k in ['ZCORN', 'COORD', 'MAPAXES', 'ACTNUM']:
if (k in attrs) and (k in sections):
if k == 'ACTNUM':
setattr(self, 'ACTNUM', sections['ACTNUM'].astype(bool))
else:
setattr(self, k, sections[k])
self.state.binary_attributes.append(k)
def _read_buffer(self, buffer, attr, logger=None, **kwargs):
if attr == 'DIMENS':
super()._read_buffer(buffer, attr, dtype=int, logger=logger, compressed=False)
elif attr in ['DX', 'DY', 'DZ', 'TOPS']:
super()._read_buffer(buffer, attr, dtype=float, logger=logger, compressed=True)
elif attr == 'ZCORN':
super()._read_buffer(buffer, attr, dtype=float, compressed=True)
elif attr == 'COORD':
super()._read_buffer(buffer, attr, dtype=float, compressed=False)
elif attr == 'MINPV':
super()._read_buffer(buffer, attr, dtype=float, compressed=False)
if 'ACTNUM' not in self.state.binary_attributes:
self._apply_minpv()
else:
if logger:
logger.info('ACTNUM is loaded from binary file: MINPV was not applied.')
elif attr in 'ACTNUM':
super()._read_buffer(buffer, attr, dtype=lambda x: bool(int(x)), logger=logger, compressed=True)
else:
super()._read_buffer(buffer, attr, logger=logger, **kwargs)
if attr in ['DX', 'DY', 'DZ']:
if np.unique(getattr(self, attr)).size != 1:
raise ValueError("Grid is not uniform ('{}').".format(attr))
setattr(self, attr, getattr(self, attr)[0])
elif attr == 'TOPS':
if np.unique(np.diff(self.TOPS)).size > 2:
raise ValueError("Grid is not uniform ('{}').".format(attr))
setattr(self, attr, self.TOPS[0])
def _apply_minpv(self):
assert hasattr(self, 'actnum')
minpv_value = self.minpv[0]
volumes = self.cell_volumes
poro = self.field.rock.poro
ntg = self.field.rock.ntg
mask = poro * volumes*ntg >= minpv_value
self.actnum = self.actnum * mask
[docs]
@apply_to_each_input
def pad_na(self, attr, fill_na=0., inplace=True):
"""Add dummy cells into the state vector in the positions of non-active cells if necessary.
Parameters
----------
attr: str, array-like
Attributes to be padded with non-active cells.
actnum: array-like of type bool
Vector representing a mask of active and non-active cells.
fill_na: float
Value to be used as filler.
inplace: bool
Modify сomponent inplace.
Returns
-------
output : component if inplace else padded attribute.
"""
_, __ = fill_na, inplace
return getattr(self, attr)
@apply_to_each_input
def _to_spatial(self, attr, inplace=True, **kwargs):
"""Spatial order 'F' transformations."""
_ = kwargs
data = getattr(self, attr)
if self.state.spatial:
return self if inplace else data
if isinstance(data, np.ndarray) and data.ndim == 1:
if attr == 'ACTNUM':
data = data.reshape(self.dimens, order='F')
elif attr == 'COORD':
nx, ny, nz = self.dimens
data = data.reshape(-1, 6)
data = data.reshape((nx + 1, ny + 1, 6), order='F')
elif attr == 'ZCORN':
nx, ny, nz = self.dimens
data = data.reshape((2, nx, 2, ny, 2, nz), order='F')
data = np.moveaxis(data, range(6), (3, 0, 4, 1, 5, 2))
data = data.reshape((nx, ny, nz, 8), order='F')
if inplace:
setattr(self, attr, data)
return self
return data
@apply_to_each_input
def _ravel(self, attr, inplace=True, **kwargs):
"""Ravel order 'F' transformations."""
_ = kwargs
data = getattr(self, attr)
if not self.state.spatial:
return self if inplace else data
if attr == 'ACTNUM':
data = data.ravel(order='F')
elif attr == 'COORD':
data = data.reshape((-1, 6), order='F').ravel()
elif attr == 'ZCORN':
nx, ny, nz = self.dimens
data = data.reshape((nx, ny, nz, 2, 2, 2), order='F')
data = np.moveaxis(data, (3, 0, 4, 1, 5, 2), range(6)).ravel(order='F')
else:
data = super()._ravel(attr=attr, order='F', inplace=False)
if inplace:
setattr(self, attr, data)
return self
return data
def _make_data_dump(self, attr, fmt=None, float_dtype=None, **kwargs):
"""Prepare data for dump."""
if fmt.upper() != 'HDF5':
return super()._make_data_dump(attr, fmt=fmt, **kwargs)
data = self.ravel(attr=attr, inplace=False)
if attr == 'ACTNUM':
return data.astype(bool)
if attr in ['ZCORN', 'COORD', 'TOPS', 'MAPAXES']:
return data if float_dtype is None else data.astype(float_dtype)
if attr == 'DIMENS':
return data.astype(int)
return data
[docs]
@state_check(lambda state: state.spatial)
def minimal_active_slices(self):
"""Get minimal cube slice that contains all active cells."""
if 'ACTNUM' in self:
pos = np.where(self.actnum)
return tuple(slice(p.min(), p.max() + 1) for p in pos)
return slice(None), slice(None), slice(None)
[docs]
@state_check(lambda state: state.spatial)
def crop_minimal_cube(self):
"""Crop to minimal cube containing active cells."""
raise NotImplementedError()
[docs]
def get_neighbors_matrix(self, connectivity=1, fill_value=-1, ravel_index=False):
"""Get indices of neighbors for all active cells.
Parameters
----------
connectivity : int, optional
Maximum number of orthogonal hops to consider a cell is
as a neighbor, by default 1.
fill_value : int, optional
Value to fill indices of inactive or absent neighbors, by default -1.
ravel_index : bool, optional
Indices in raveled form, by default False.
Returns
-------
res : misc
Matrix of active neighbors and matrix of distances if 'calculate_distances'.
"""
actnum = self.to_spatial(attr='ACTNUM', inplace=False)
res, invalid_cells = get_connectivity_matrix(actnum, connectivity)
if ravel_index:
res = np.ravel_multi_index((res[..., 0], res[..., 1], res[..., 2]), self.dimens, order='F')
res[invalid_cells] = fill_value
return res
[docs]
def calculate_neighbours_distances(self, connectivity=1, fill_value=-1, neighbours_matrix=None):
"""Calculate distances between neighbors for all active cells.
Parameters
----------
connectivity : int, optional
Maximum number of orthogonal hops to consider a cell is
as a neighbor., by default 1.
fill_value : int, optional
Value to fill indices of inactive or absent neighbors, by default -1.
neighbours_matrix: numpy.ndarray, optional
Matrix with indices of neighbors in unraveled form.
Returns
-------
numpy.ndarray
Matrix of distances.
"""
actnum = self.to_spatial(attr='ACTNUM', inplace=False)
if neighbours_matrix is None:
neighbours_matrix = self.get_neighbors_matrix(
connectivity=connectivity,
fill_value=fill_value,
ravel_index=False
)
invalid_cells = np.where(
np.any(neighbours_matrix == fill_value, axis=-1),
np.ones(shape=neighbours_matrix.shape[:2]),
np.zeros(shape=neighbours_matrix.shape[:2])
).astype(bool)
neighbor_centers = self.cell_centroids[
neighbours_matrix[..., 0], neighbours_matrix[..., 1], neighbours_matrix[..., 2]
]
active_cells_centers = self.cell_centroids[actnum]
distances = np.linalg.norm(
neighbor_centers - np.tile(active_cells_centers[:, np.newaxis, :],
(1, neighbor_centers.shape[1], 1)),
axis=2
)
distances[invalid_cells] = fill_value
return distances
[docs]
class CornerPointGrid(Grid):
"""Corner point grid."""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
if 'MAPAXES' not in self:
setattr(self, 'MAPAXES', np.array([0, 1, 0, 0, 1, 0]))
self.to_spatial()
@property
def origin(self):
"""Grid axes origin relative to the map coordinates."""
return np.array([self.mapaxes[2], self.mapaxes[3], self.zcorn[0, 0, 0, 0]])
@cached_property(
lambda self, x: x.reshape(tuple(self.dimens) + (8, 3), order='F') if self.state.spatial
else x.reshape(-1, 8, 3, order='F'), modify_cache=True
)
def xyz(self):
"""x, y, z coordinates of cells vertices."""
if not self.state.spatial:
zcorn = self._to_spatial(attr='ZCORN', inplace=False)
coord = self._to_spatial(attr='COORD', inplace=False)
return numba_get_xyz(self.dimens, zcorn, coord).reshape([-1, 8, 3], order='F')
return numba_get_xyz(self.dimens, self.zcorn, self.coord)
@cached_property(
lambda self, x: x.reshape(tuple(self.dimens) + (3,), order='F') if self.state.spatial
else x.reshape(-1, 3, order='F'), modify_cache=True
)
def xyz_min(self):
"""Minimums of x, y, z coordinates of cells vertices."""
return np.min(self.xyz, axis=-2)
@cached_property(
lambda self, x: x.reshape(tuple(self.dimens) + (3,), order='F') if self.state.spatial
else x.reshape(-1, 3, order='F'), modify_cache=True
)
def xyz_max(self):
"""Maximums of x, y, z coordinates of cells vertices."""
return np.max(self.xyz, axis=-2)
@cached_property(
lambda self, x: x.reshape(tuple(self.dimens) + (3,), order='F') if self.state.spatial
else x.reshape(-1, 3, order='F'), modify_cache=True
)
def cell_centroids(self):
"""Centroids of cells."""
return self.xyz.mean(axis=-2)
@cached_property
@state_check(lambda state: state.spatial)
def cell_volumes(self):
"""Volumes of cells."""
return numba_get_volumes(self.xyz)
@cached_property
@state_check(lambda state: state.spatial)
def bounding_box(self):
"""Pair of diagonal corner points for grid's bounding box."""
min_vals = self.xyz_min.min(axis=(0, 1, 2))
max_vals = self.xyz_max.max(axis=(0, 1, 2))
return np.array([min_vals, max_vals])
[docs]
@state_check(lambda state: state.spatial)
def minimal_active_bounds(self):
"""Get z coordinates of top and bottom bounds of active cells."""
zcorn = self.zcorn.reshape(self.zcorn.shape[:3] + (2, 4))
z_top = get_top_z_coords(zcorn, self.actnum)
z_bottom = -get_top_z_coords(-zcorn[:, :, ::-1, ::-1], self.actnum[:, :, ::-1])
return z_top, z_bottom
[docs]
@state_check(lambda state: state.spatial)
def upscale(self, factors=(2, 2, 2), actnum_upscale='vote'):
"""Upscale grid according to factors given.
Parameters
----------
factors : tuple, int
Scale factors along each axis. If int, factors are the same for each axis.
actnum_upscale : str
Method to actnum upscaling. If 'vote', upscaled cell is active if majority
of finer cells are active. If 'any', upscaled cell is active if any
of finer cells is active. Default to 'vote'.
Returns
-------
grid : CornerPointGrid
Upscaled grid.
"""
factors = np.atleast_1d(factors)
if factors.size == 1:
factors = np.repeat(factors, 3)
coord = self.coord[::factors[0], ::factors[1]]
dimens = self.dimens // factors
d0, d1, d2 = dimens * factors
s0, s1, s2 = factors
zcorn = np.zeros(tuple(dimens) + (8,), dtype=self.zcorn.dtype)
zcorn[:, :, :, 0] = self.zcorn[:d0:s0, :d1:s1, :d2:s2, 0]
zcorn[:, :, :, 1] = self.zcorn[s0 - 1:d0:s0, :d1:s1, :d2:s2, 1]
zcorn[:, :, :, 2] = self.zcorn[:d0:s0, s1 - 1:d1:s1, :d2:s2, 2]
zcorn[:, :, :, 3] = self.zcorn[s0 - 1:d0:s0, s1 - 1:d1:s1, :d2:s2, 3]
zcorn[:, :, :, 4] = self.zcorn[:d0:s0, :d1:s1, s2 - 1:d2:s2, 4]
zcorn[:, :, :, 5] = self.zcorn[s0 - 1:d0:s0, :d1:s1, s2 - 1:d2:s2, 5]
zcorn[:, :, :, 6] = self.zcorn[:d0:s0, s1 - 1:d1:s1, s2 - 1:d2:s2, 6]
zcorn[:, :, :, 7] = self.zcorn[s0 - 1:d0:s0, s1 - 1:d1:s1, s2 - 1:d2:s2, 7]
if 'ACTNUM' not in self:
grid = self.__class__(coord=coord, dimens=dimens, zcorn=zcorn, mapaxes=self.mapaxes)
grid.set_state(spatial=True)
return grid
out = rolling_window(self.actnum, factors)
if actnum_upscale == 'vote':
actnum = np.mean(out, axis=(-3, -2, -1)) > 0.5
elif actnum_upscale == 'any':
actnum = np.sum(out, axis=(-3, -2, -1)) > 0
else:
raise ValueError('Unknown mode of actnum upscaling: {}.'.format(actnum_upscale))
grid = self.__class__(coord=coord, dimens=dimens, zcorn=zcorn,
mapaxes=self.mapaxes, actnum=actnum)
grid.set_state(spatial=True)
return grid
[docs]
@state_check(lambda state: state.spatial)
def orthogonalize(self, dimens, only_active=False):
"""Construct orthogonal uniform grid. Actnum tranfer should be made separately."""
nx, ny, nz = dimens
xyz_corn = self.xyz[self.actnum] if only_active else self.xyz
x_min = xyz_corn[..., 0].min()
x_max = xyz_corn[..., 0].max()
y_min = xyz_corn[..., 1].min()
y_max = xyz_corn[..., 1].max()
z_min = xyz_corn[..., 2].min()
z_max = xyz_corn[..., 2].max()
dx = (x_max - x_min) / nx
dy = (y_max - y_min) / ny
dz = (z_max - z_min) / nz
mapaxes = np.array([x_min, y_min + 100, x_min, y_min, x_min + 100, y_min])
grid = OrthogonalUniformGrid(dimens=dimens, dx=dx, dy=dy, dz=dz,
tops=z_min, mapaxes=mapaxes)
grid.set_state(spatial=True)
return grid
[docs]
@state_check(lambda state: state.spatial)
def crop_minimal_cube(self):
"""Crop to minimal cube containing active cells.
Returns
-------
(grid, min_slices) : tuple
New grid and slices for active cube.
"""
min_slices = self.minimal_active_slices()
actnum = self.actnum[min_slices]
zcorn = self.zcorn[min_slices]
min_coord_slices = tuple(slice(s.start, s.stop + 1) for s in min_slices[:2])
coord = self.coord[min_coord_slices]
grid = self.__class__(actnum=actnum, coord=coord,
dimens=np.array(actnum.shape),
zcorn=zcorn, mapaxes=self.mapaxes)
grid.set_state(spatial=True)
return grid, min_slices
[docs]
@state_check(lambda state: state.spatial)
def crop_minimal_grid(self, nz, fillna=None):
"""Create a new grid in a region between upper and bottom surfaces enclosing active cells.
Actnum tranfer should be made separately.
Parameters
----------
nz : int
Third dimension of the cunstructed grid.
fillna : scalar
Filling value for nan coordinates.
Returns
-------
(grid, grid_mask, z_top, z_bottom) : tuple
"""
z_top, z_bottom = self.minimal_active_bounds()
with warnings.catch_warnings(): # ignore comparison with np.nan
warnings.simplefilter("ignore")
grid_mask = ((self.zcorn[:, :, :, 0] >= np.expand_dims(z_top[:-1, :-1], -1)) &
(self.zcorn[:, :, :, 4] <= np.expand_dims(z_bottom[:-1, :-1], -1)))
def _fillna(x):
"""Replace np.nan with value if value is given."""
if not np.isnan(x):
return x
return x if fillna is None else fillna
grid_points = np.stack([np.linspace(_fillna(start), _fillna(stop), nz + 1)
for start, stop in zip(z_top.ravel(), z_bottom.ravel())])
grid_points = grid_points.reshape(z_top.shape + (nz + 1,))
zcorn = np.stack([grid_points[:-1, :-1, :-1], # 0
grid_points[1:, :-1, :-1], # 1
grid_points[:-1, 1:, :-1], # 2
grid_points[1:, 1:, :-1], # 3
grid_points[:-1, :-1, 1:], # 4
grid_points[1:, :-1, 1:], # 5
grid_points[:-1, 1:, 1:], # 6
grid_points[1:, 1:, 1:]], axis=-1)
grid = self.__class__(coord=self.coord, dimens=np.array(zcorn.shape[:3]),
zcorn=zcorn, mapaxes=self.mapaxes)
grid.set_state(spatial=True)
return grid, grid_mask, z_top, z_bottom
@property
def as_corner_point(self):
"""Returns itself."""
return self
@state_check(lambda state: state.spatial)
def create_vtk_grid(self, use_only_active=True, scaling=True, **kwargs):
"""Creates pyvista unstructured grid object.
Returns
-------
grid : `pyvista.core.pointset.UnstructuredGrid` object
"""
_ = kwargs
self._vtk_grid_params = {'use_only_active': use_only_active, 'cell_size': None, 'scaling': scaling}
cells = self.xyz
indexes = np.moveaxis(np.indices(self.dimens), 0, -1)
if 'ACTNUM' in self and use_only_active:
cells = cells[self.actnum]
indexes = indexes[self.actnum]
else:
cells = cells.reshape((-1,) + cells.shape[3:])
indexes = indexes.reshape((-1,) + indexes.shape[3:])
self._cell_id_d = dict(enumerate(indexes))
cells[:, [2, 3]] = cells[:, [3, 2]]
cells[:, [6, 7]] = cells[:, [7, 6]]
n_cells = cells.shape[0]
cells = cells.reshape(-1, cells.shape[-1], order='C')
cells = numpy_support.numpy_to_vtk(cells, deep=True)
points = vtk.vtkPoints()
points.SetNumberOfPoints(8*n_cells)
points.SetData(cells)
cell_array = vtk.vtkCellArray()
cell = vtk.vtkHexahedron()
connectivity = np.insert(range(8 * n_cells), range(0, 8 * n_cells, 8), 8).astype(np.int64)
cell_array.SetCells(n_cells, numpy_support.numpy_to_vtkIdTypeArray(connectivity, deep=True))
self._vtk_grid = vtk.vtkUnstructuredGrid()
self._vtk_grid.SetPoints(points)
self._vtk_grid.SetCells(cell.GetCellType(), cell_array)
@state_check(lambda state: state.spatial)
def create_vtk_locator(self, use_only_active=True, scaling=False, **kwargs):
"""Creates locator and mapping dictionary.
Returns
-------
grid : `pyvista.core.pointset.UnstructuredGrid` object
"""
_ = kwargs
grid_params = {'use_only_active': use_only_active, 'cell_size': None, 'scaling': scaling}
if self._vtk_grid is None or self._vtk_grid_params != grid_params:
self.create_vtk_grid(**grid_params)
geo_filter = vtk.vtkGeometryFilter()
geo_filter.SetInputData(self._vtk_grid)
geo_filter.Update()
polydata = geo_filter.GetOutput()
locator = vtk.vtkModifiedBSPTree()
# locator.LazyEvaluationOff() #Deprecated since vtk==9.2.0
locator.SetDataSet(polydata)
locator.AutomaticOn()
locator.BuildLocator()
self._vtk_locator = locator
@state_check(lambda state: state.spatial)
def point_inside_cell(self, point, cell_idx, tolerance=1e-8):
"""Determines whether point is inside cell or not.
Returns
-------
inside : bool
Notes
-----
Result might be inconsistent due to the stochastic nature of the algorithm.
Decreasing the tolerance might help.
"""
x = self.xyz[cell_idx[0], cell_idx[1], cell_idx[2]].copy()
x[[2, 3]] = x[[3, 2]]
x[[6, 7]] = x[[7, 6]]
pts = [(0, 1, 2, 3), (4, 5, 6, 7), (0, 1, 5, 4),
(1, 2, 6, 5), (2, 3, 7, 6), (3, 0, 4, 7)]
cube = vtk.vtkPolyData()
points = vtk.vtkPoints()
points.SetDataTypeToDouble()
polys = vtk.vtkCellArray()
for i in range(8):
points.InsertPoint(i, x[i])
for i in range(6):
polys.InsertNextCell(mk_vtk_id_list(pts[i]))
cube.SetPoints(points)
cube.SetPolys(polys)
points = np.array([point])
vtk_pts = vtk.vtkPoints()
vtk_pts.SetDataTypeToDouble()
vtk_pts.SetData(numpy_support.numpy_to_vtk(points, deep=1))
# Make the poly data
poly_pts_vtp = vtk.vtkPolyData()
poly_pts_vtp.SetPoints(vtk_pts)
enclosed_points_filter = vtk.vtkSelectEnclosedPoints()
enclosed_points_filter.SetTolerance(tolerance)
enclosed_points_filter.SetSurfaceData(cube)
enclosed_points_filter.SetInputData(poly_pts_vtp)
enclosed_points_filter.Update()
return enclosed_points_filter.IsInside(0)
[docs]
def faces_centers(self, cells_indices):
"""Calculate coordinates of cell faces centers.
Parameters
----------
cells_indices : List[np.ndarray]
Indices of the cells.
Returns
-------
np.ndarray
coordinates of cell faces centers.
"""
xyz = self.xyz[cells_indices[0], cells_indices[1], cells_indices[2]]
centers = []
for vertices in (
(0, 2, 4, 6), # left
(1, 3, 5, 7), # right
(0, 1, 4, 5), # back
(2, 3, 6, 7), # front
(0, 1, 2, 3), # top
(4, 5, 6, 7) # bottom
):
centers.append(xyz[..., vertices, :].mean(axis=-2))
centers = np.array(centers)
if centers.ndim == 3:
return np.moveaxis(centers, (0, 1), (1, 0))
return centers
[docs]
def cell_sizes(self, cell_indices):
"""Calculate approximate sizes of cells.
Parameters
----------
cells_indices : List[np.ndarray]
Indices of the cells.
Returns
-------
np.ndarray
Sizes of cells.
"""
faces_centers = self.faces_centers(cell_indices)
n_cells = cell_indices[0].size if isinstance(cell_indices[0], np.ndarray) else 1
sizes = np.zeros((n_cells, 3))
for i in range(3):
sizes[:, i] = np.sqrt((((faces_centers[..., 2*i+1, :] - faces_centers[..., 2*i+0, :])**2)).sum(axis=-1))
return sizes
[docs]
def cell_bases(self, cell_indices):
"""Calculate basis vectors of coordinate systems connected to cells.
Parameters
----------
cells_indices : List[np.ndarray]
Indices of the cells.
Returns
-------
np.ndarray
Basis vectors of coordinate systems connected to cells.
"""
faces_centers = self.faces_centers(cell_indices)
i = faces_centers[..., 1, :] - faces_centers[..., 0, :]
j = faces_centers[..., 3, :] - faces_centers[..., 2, :]
z = np.cross(i, j)
hat3 = (z.T / np.linalg.norm(z, axis=-1)).T
x = i + np.cross(j, hat3)
y = j - np.cross(i, hat3)
hat1 = (x.T / np.linalg.norm(x, axis=-1)).T
hat2 = (y.T / np.linalg.norm(y, axis=-1)).T
return np.moveaxis(np.stack((hat1, hat2, hat3)), range(2), (1, 0))
[docs]
@state_check(lambda state: state.spatial)
def map_grid(self):
"""Map pillars (`COORD`) to axis defined by `MAPAXES'.
Returns
-------
CornerPointGrid
Grid with updated `COORD` and `MAPAXES` fields.
"""
if not np.isclose(self.ex.dot(self.ey), 0):
raise ValueError('`ex` and `ey` vectors should be orthogonal.')
new_basis = np.vstack((self.ex, self.ey)).T
self.coord[:, :, :2] = self.coord[:, :, :2].dot(new_basis) + self.origin[:2]
self.coord[:, :, 3:5] = self.coord[:, :, 3:5].dot(new_basis) + self.origin[:2]
setattr(self, 'MAPAXES', np.array([0, 1, 0, 0, 1, 0]))
return self
def specify_grid(grid):
"""Specify grid class: `CornerPointGrid` or `OrthogonalUniformGrid`.
Parameters
----------
grid : Grid
Initial grid.
Returns
-------
CornerPointGrid or OrthogonalUniformGrid
specified grid.
"""
if not isinstance(grid, (CornerPointGrid, OrthogonalUniformGrid)):
if ('DX' in grid) and ('DY' in grid) and ('DZ' in grid):
grid = OrthogonalUniformGrid(**dict(grid.items()), field=grid.field)
else:
grid = CornerPointGrid(**dict(grid.items()), field=grid.field)
return grid