"""Grid class for the pygeon package."""
from typing import Tuple
import numpy as np
import porepy as pp
import scipy.sparse as sps
import pygeon as pg
"""
Acknowledgments:
The functionalities related to the ridge computations are modified from
github.com/anabudisa/md_aux_precond developed by Ana Budiša and Wietse M. Boon.
"""
[docs]
class Grid(pp.Grid):
"""
Grid class represents a geometric grid object, in addition to the pp.Grid class it
implements the following attributes and methods.
"""
[docs]
def __init__(self, *args, **kwargs) -> None:
"""
Initialize a Grid object.
Args:
*args: Variable length argument list.
**kwargs: Arbitrary keyword arguments.
Returns:
None
"""
super().__init__(*args, **kwargs)
self.face_nodes: sps.csc_array
self.cell_faces: sps.csc_array
[docs]
def compute_geometry(self) -> None:
"""
Defines grid entities of codim 2 and 3.
The entities are referred to by their codimension:
0: "cells"
1: "faces"
2: "ridges"
3: "peaks"
This method computes the geometry of the grid by calling the
superclass's compute_geometry method, computing the ridge
and peak connectivities, and storing the edge lengths and mesh size.
Args:
None
Returns:
None
"""
super().compute_geometry()
self.compute_rotation_matrix()
self.compute_ridges()
self.compute_edge_properties()
self.compute_mesh_size()
[docs]
def compute_ridges(self) -> None:
"""
Computes the ridges of the grid and assigns the following attributes:
- num_ridges: number of ridges
- num_peaks: number of peaks
- face_ridges: connectivity between each face and ridge
- ridge_peaks: connectivity between each ridge and peak
- tags['tip_ridges']: tags for entities at fracture tips
- tags['tip_peaks']: tags for entities at fracture tips
Args:
None
Returns:
None
"""
if self.dim == 3:
self._compute_ridges_3d()
elif self.dim == 2:
self._compute_ridges_2d()
else: # The grid is of dimension 0 or 1.
self._compute_ridges_01d()
self.tag_ridges()
def _compute_ridges_01d(self) -> None:
"""
Assign the number of ridges, number of peaks, and connectivity matrices to a
grid of dimension 0 or 1.
This method calculates the number of ridges and peaks in a grid of dimension
0 or 1.
It also initializes the ridge_peaks and face_ridges matrices with the
appropriate dimensions.
Args:
None
Returns:
None
"""
self.num_peaks = 0
self.num_ridges = 0
self.ridge_peaks = sps.csc_array((self.num_peaks, self.num_ridges), dtype=int)
self.face_ridges = sps.csc_array((self.num_ridges, self.num_faces), dtype=int)
def _compute_ridges_2d(self) -> None:
"""
Assign the number of ridges, number of peaks, and connectivity matrices to a
grid of dimension 2.
This method computes the number of ridges, number of peaks, and connectivity
matrices for a 2-dimensional grid. It also computes the face-ridge orientation
based on the rotated normal and the difference vector between the ridges.
Args:
None
Returns:
None
"""
self.num_peaks = 0
self.num_ridges = self.num_nodes
self.ridge_peaks = sps.csc_array((self.num_peaks, self.num_ridges), dtype=int)
# We compute the face tangential by mapping the face normal to a reference grid
# in the xy-plane, rotating locally, and mapping back.
R = self.rotation_matrix
loc_rot = np.array([[0.0, -1.0], [1.0, 0.0]])
rot = R.T @ loc_rot @ R
rotated_normal = rot @ self.face_normals
# The face-ridge orientation is determined by whether the rotated normal
# coincides with the difference vector between the ridges.
face_ridges = self.face_nodes.copy().astype(int)
face_ridges.data[::2] *= -1
face_tangents = self.nodes @ face_ridges
orients = np.sign(np.sum(rotated_normal * face_tangents, axis=0))
self.face_ridges = face_ridges @ sps.diags_array(orients)
def _compute_ridges_3d(self) -> None:
"""
Assign the number of ridges, number of peaks, and connectivity matrices to a
grid of dimension 3.
This method computes the number of ridges, number of peaks, and connectivity
matrices for a 3-dimensional grid. It calculates the ridges between each pair of
nodes in each face, determines the orientation of each ridge with respect to the
face, and generates the ridge-peak and face-ridge connectivity matrices.
Args:
None
Returns:
None
"""
self.num_peaks = self.num_nodes
# Pre-allocation
ridges = np.empty((2, self.face_nodes.nnz), dtype=int)
# Define ridges between each pair of nodes assuming ordering in face_nodes is
# done according to right-hand rule
ridges[0] = self.face_nodes.indices
ridges[1] = np.roll(ridges[0], -1)
ridges[1, self.face_nodes.indptr[1:] - 1] = ridges[
0, self.face_nodes.indptr[:-1]
]
# There are as many face-ridge pairs as face-node pairs
fr_indptr = self.face_nodes.indptr
# Save orientation of each ridge w.r.t. the face
orientations = np.sign(ridges[1, :] - ridges[0, :])
# Ridges are oriented from low to high node indices, i.e. [0,1], [0,2], [1,2].
ridges.sort(axis=0)
# Identify the ridges based on unique pairs of peaks. We do this by creating a
# sparse array. If ridge r has nodes (x_i, x_j) then A_ij = index_r. This way we
# can easily look up the unique index of a ridge based on its nodes.
ridge_index = sps.coo_array(
(np.ones(ridges.shape[1]), (ridges[0], ridges[1])), dtype=int
)
ridge_index.sum_duplicates()
ridge_index.data = np.arange(ridge_index.nnz)
# Extract the indices of the found ridges
indices = ridge_index.tocsr()[ridges[0], ridges[1]]
# We can now replace the ridges array to the one without duplicates
ridges = np.vstack((ridge_index.row, ridge_index.col))
self.num_ridges = np.size(ridges, 1)
# Generate ridge-peak connectivity such that
# ridge_peaks(i, j) = +/- 1:
# ridge j points to/away from peak i
indptr = np.arange(0, ridges.size + 1, 2)
ind = np.ravel(ridges, order="F")
data = np.ones(ridges.size)
data[::2] *= -1
self.ridge_peaks = sps.csc_array((data, ind, indptr), dtype=int)
# Generate face_ridges such that
# face_ridges(i, j) = +/- 1:
# face j has ridge i with same/opposite orientation
# with the orientation defined according to the right-hand rule
self.face_ridges = sps.csc_array((orientations, indices, fr_indptr), dtype=int)
[docs]
def tag_ridges(self) -> None:
"""
Tag the peaks and ridges of the grid located on fracture tips.
This method tags the peaks and ridges of the grid that are located on fracture
tips. It sets the "tip_peaks" and "tip_ridges" tags in the grid object.
For 2D grids, the "tip_ridges" tag is determined based on the "tip_faces" tag
and the face ridges.
For 3D grids, the "tip_ridges" tag is initialized as an array of zeros.
The "domain_boundary_ridges" tag is also set based on the face ridges and the
"domain_boundary_faces" tag.
Args:
None
Returns:
None
"""
self.tags["tip_peaks"] = np.zeros(self.num_peaks, dtype=bool)
fr_bool = self.face_ridges.astype("bool")
if self.dim == 2:
self.tags["tip_ridges"] = fr_bool @ self.tags["tip_faces"]
else:
self.tags["tip_ridges"] = np.zeros(self.num_ridges, dtype=bool)
bd_ridges = fr_bool @ self.tags["domain_boundary_faces"]
self.tags["domain_boundary_ridges"] = bd_ridges.astype(bool)
[docs]
def compute_subvolumes(
self, return_subsimplices: bool = False
) -> Tuple[sps.csc_array, sps.csc_array] | sps.csc_array:
"""
Compute the subvolumes of the grid.
Args:
return_subsimplices (bool, optional): Whether to return the sub-simplices.
Defaults to False.
Returns:
sps.csc_array: The computed subvolumes with each entry [node, cell]
describing the signed measure of the associated sub-volume
"""
sub_simplices = sps.csc_array(self.cell_faces.copy().astype(float))
faces, cells, orient = sps.find(self.cell_faces)
normals = self.face_normals[:, faces] * orient
rays = self.face_centers[:, faces] - self.cell_centers[:, cells]
sub_simplices[faces, cells] = np.sum(normals * rays, axis=0) / self.dim
nodes_per_face = np.array(self.face_nodes.sum(axis=0)).flatten()
div_by_nodes_per_face = sps.diags_array(1.0 / nodes_per_face)
subv: sps.csc_array = self.face_nodes @ div_by_nodes_per_face @ sub_simplices
if return_subsimplices:
return subv, sub_simplices
else:
return subv
[docs]
def compute_opposite_nodes(self, recompute=False) -> sps.csc_array:
"""
Computes a matrix containing the index of the opposite node
for every (face, cell) pair. Sets it as an attribute for later use.
Args:
recompute (bool, optional): Whether to recompute the opposite nodes.
Defaults to False.
Returns:
sps.csc_array: The index k of the opposite node is in the entry (face, cell)
"""
if recompute or not hasattr(self, "opposite_nodes"):
cell_nodes = self.cell_nodes()
if self.dim == 0:
self.opposite_nodes = sps.csc_array((0, 1), dtype=int)
return self.opposite_nodes
if not np.all(cell_nodes.sum(axis=0) == self.dim + 1):
raise NotImplementedError(
"Grid is not simplicial; cannot compute opposite node."
)
faces, cells, _ = sps.find(self.cell_faces)
opposites = cell_nodes[:, cells] - self.face_nodes[:, faces].astype(bool)
self.opposite_nodes = sps.csc_array((opposites.indices, (faces, cells)))
return self.opposite_nodes
[docs]
def compute_edge_properties(self) -> None:
"""
Computes and stores the tangent vectors and lengths
of the grid edges (1D entities).
Args:
None
Returns:
None
"""
match self.dim:
case 0:
self.edge_tangents = np.zeros((0, pg.AMBIENT_DIM))
self.edge_lengths = np.zeros(0)
self.num_edges = 0
return
case 1:
edge_nodes = self.cell_faces
case 2:
edge_nodes = self.face_ridges
case 3:
edge_nodes = self.ridge_peaks
self.edge_tangents = self.nodes @ edge_nodes
self.edge_lengths = np.sqrt(np.sum(self.edge_tangents**2, axis=0))
self.num_edges = self.edge_lengths.size
[docs]
def compute_mesh_size(self) -> None:
"""
Computes and stores the typical
mesh size as the mean of the edge lengths.
Args:
None
Returns:
None
"""
if self.dim == 0:
self.mesh_size = 0.0
else:
self.mesh_size = float(np.mean(self.edge_lengths))
[docs]
def compute_rotation_matrix(self) -> None:
"""
Computes and stores the rotation matrix that maps the subdomain to the xy-plane
or x-axis.
Args:
None
Returns:
None
"""
if np.any(self.nodes[self.dim :, :]):
*_, R, keep_dims, _ = pp.map_geometry.map_grid(self)
self.rotation_matrix = R[keep_dims, :]
else:
self.rotation_matrix = np.eye(self.dim, pg.AMBIENT_DIM)
[docs]
def copy(self):
"""Create a new instance with some attributes deep-copied from the grid.
Returns:
A deep copy of ``self``. Some predefined attributes are also copied.
"""
h = super().copy()
pg.convert_from_pp(h)
return h