pygeon.discretizations.fvm.tpsa module
Module for two-point stress approximation discretization.
- class pygeon.discretizations.fvm.tpsa.TPSA(keyword='unitary_data')[source]
Bases:
FiniteVolumeDiscretizationA vectorized implementation of the two-point stress approximation method for elasticity of Nordbotten and Keilegavlen (2025).
Degrees of freedom are given by the cell center values and the variables are ordered as [ux, uy, uz, rx, ry, rz, p], with the cell index varying fastest.
Our implementation differs from Porepy (v1.13) in the imposition of boundary conditions. In particular, we made the following changes:
The order of R and Xi_tilde in the [0, 1] block of (A.2.25)
The order of n and Xi_tilde in the [0, 2] block of (A.2.25)
Changed “delta R^2” to “R delta R” in the [1, 1] block of (A.2.25)
Used the delta^mu_k in the normal direction in the [2, 2] block of (A.2.25)
These adaptations are necessary for consistency with rolling boundary conditions. We moreover used signed distances between face and cell centers.
Equation numbers in the comments and docstrings refer to the manuscript: https://doi.org/10.1016/j.camwa.2025.07.035
- __init__(keyword='unitary_data')[source]
Initialize the TPSA object.
- Parameters:
keyword (str) – The keyword used to identify the discretization method. Default is pg.UNITARY_DATA.
- Returns:
None
- ndof_per_cell(sd)[source]
Returns the number of degrees of freedom per cell.
- Parameters:
sd (pg.Grid) – The grid object.
- Returns:
The number of degrees of freedom per cell.
- Return type:
- interpolate(sd, displacement, rotation, solid_pressure)[source]
Interpolate a triplet of functions onto the finite volume space.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
displacement (Callable) – A function that returns the displacement values at coordinates.
rotation (Callable) – A function that returns the rotation values at coordinates.
solid_pressure (Callable) – A function that returns the solid pressure values at coordinates.
- Returns:
The values of the degrees of freedom
- Return type:
np.ndarray
- assemble_accumulation_terms(sd, data)[source]
Assemble the zeroth-order terms on the diagonal of (3.9).
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
lame_mu (np.ndarray) – The Lamé parameter mu
lame_lambda (np.ndarray) – The Lamé parameter lambda
- Returns:
The diagonal mass matrix
- Return type:
sps.csc_array
- precompute_arrays(sd, data=None)[source]
Precomputations on the grid for easy access later.
This function is typically called twice, once for the left-hand side, and once for the right.
- compute_weighted_dists(sd, data, find_cell_faces)[source]
Computes delta_k^i / mu_i from (2.1) for every physical face-cell pair (k, i). Boundary conditions are handled later.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
data (dict) – The data dictionary.
find_cell_faces (Tuple) – Output of scipy.sparse.find on sd.cell_faces.
- Returns:
The weighted distances.
- Return type:
np.ndarray
- extend_faces_and_distances(sd, data, faces, weighted_dists)[source]
Incorporate the boundary conditions by extending the face and distance arrays.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
data (dict) – The data dictionary.
faces (np.ndarray) – The array of face indices
weighted_dists (np.ndarray) – The array of weighted distances
- Returns:
The extended array of face indices. np.ndarray: The extended array of weighted distances.
- Return type:
np.ndarray
- compute_delta_mu_k(faces, dists)[source]
Compute the delta^mu_k of (3.5) given by 0.5 * ( mu_i delta_k^-i + mu_j delta_k^-j)^-1 for each face k with neighboring cells (i,j).
- Parameters:
faces (np.ndarray) – The extended array of faces.
dists (np.ndarray) – The extended array of weighted distances.
- Returns:
The array of $delta_k^mu$.
- Return type:
np.ndarray
- compute_harmonic_avg(faces, dists)[source]
Compute the harmonic average of mu from (3.5), divided by delta_k, at each face: mu_effective = ( delta_k^i / mu_i + delta_k^j / mu_j)^-1
- Parameters:
faces (np.ndarray) – The extended array of faces.
dists (np.ndarray) – The extended array of weighted distances.
- Returns:
The face-wise harmonic average of $mu$.
- Return type:
np.ndarray
- assemble_dual_var_map(sd, data)[source]
Assemble the mapping from cell-based primary variables to face-based dual variables.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
- Returns:
The matrix mapping primary to dual variables
- Return type:
sps.csc_array
- assemble_rot(sd)[source]
The operator R^n that performs a cross product with the normal vector n.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
- Returns:
The R^n matrix
- Return type:
sps.csc_array
- assemble_ndot(sd)[source]
The operator that performs a dot product with the normal vector n.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
- Returns:
The n cdot matrix
- Return type:
sps.csc_array
- assemble_rot_rot_bdry_terms(sd, cached_arrays)[source]
The operator R^n delta R^n that is on the [1, 1] block of (A2.25).
There is a slight discrepancy with the paper, because a simpler class of boundary conditions are assumed there. This implementation is the generalization to more involved boundary conditions, such as rollers.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
cached_arrays (dict) – The output of self.precompute_arrays
- Returns:
The double rotation matrix
- Return type:
sps.csc_array
- assemble_Xi(cached_arrays)[source]
Compute the averaging operator Xi from (2.5).
Displacement bc are handled by delta_mu_k = 0. Traction bc are handled since 2 * delta_mu_k * mu / delta = 1. Spring bc are handled because the spring constant is contained in delta_mu_k.
- convert_to_xi_tilde_inplace(Xi)[source]
Compute the complementary averaging operator Xi_tilde from (2.6). NOTE: This is an in-place operation.
- assemble_first_column(sd, Xi_list)[source]
Assemble the off-diagonal terms in the first column of (3.7). These are computed together because their construction uses similar components.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
Xi_list (list) – The averaging operators in the coordinate directions
- Returns:
The operator that averages and then crosses with n n_xi (sps.csc_array): The operator that averages and then dots with n
- Return type:
R_xi (sps.csc_array)
- assemble_first_row(sd, Xi_list)[source]
Assemble the off-diagonal terms in the first row of (3.7). These are computed together because their construction uses similar components.
This is a generalization compared to the paper to handle more involved boundary conditions. In particular, we have to change the order of the operators
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
Xi_list (list) – The tilde averaging operators in the coordinate directions
- Returns:
The operator that crosses with n and then averages n_xi (sps.csc_array): The operator that multiplies with n and then averages
- Return type:
R_xi (sps.csc_array)
- assemble_bdry_dual_var_map(sd, data=None)[source]
Assemble the second matrix on the right-hand side of (A2.25).
Slight generalization: the [1, 0] and [2, 0] blocks first weigh with delta and then rot/dot with n.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
data (dict) – The data dictionary
- Returns:
the matrix to be multiplied with the boundary data g
- Return type:
sps.csc_array
- assemble_body_force(sd, func)[source]
Assemble the right-hand side for a given body-force function.
- Parameters:
sd (pg.Grid) – Grid, or a subclass.
func (Callable) – The body force function.
- Returns:
The right-hand side vector.
- Return type:
np.ndarray