API Reference

This section provides the API reference for the key modules and functions.

common

common.block_diag(mats, format=None, dtype=None)

Build a block diagonal sparse matrix or array from provided matrices.

Parameters:

• mats (sequence of matrices or arrays) – Input matrices or arrays.

• format (str, optional) – The sparse format of the result (e.g., “csr”). If not given, the result is returned in “coo” format.

• dtype (dtype specifier, optional) – The data-type of the output. If not given, the dtype is determined from that of blocks.

Returns:

res – If at least one input is a sparse array, the output is a sparse array. Otherwise the output is a sparse matrix.

Return type:

sparse matrix or array

Notes

Added in version 0.11.0.

See also

block_array, diags_array

Examples

>>> from scipy.sparse import coo_array, block_diag
>>> A = coo_array([[1, 2], [3, 4]])
>>> B = coo_array([[5], [6]])
>>> C = coo_array([[7]])
>>> block_diag((A, B, C)).toarray()
array([[1, 2, 0, 0],
       [3, 4, 0, 0],
       [0, 0, 5, 0],
       [0, 0, 6, 0],
       [0, 0, 0, 7]])

common.bmat(blocks, format=None, dtype=None)

Build a sparse array or matrix from sparse sub-blocks

Note: block_array is preferred over bmat. They are the same function except that bmat returns a deprecated sparse matrix when none of the inputs are sparse arrays.

Warning

This function returns a sparse matrix when no inputs are sparse arrays. You are encouraged to use block_array to take advantage of the sparse array functionality.

Parameters:

• blocks (array_like) – Grid of sparse matrices with compatible shapes. An entry of None implies an all-zero matrix.

• format ({'bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil'}, optional) – The sparse format of the result (e.g. “csr”). By default an appropriate sparse matrix format is returned. This choice is subject to change.

• dtype (dtype, optional) – The data-type of the output matrix. If not given, the dtype is determined from that of blocks.

Returns:

bmat – If any block in blocks is a sparse array, return a sparse array. Otherwise return a sparse matrix.

If you want a sparse array built from blocks that are not sparse arrays, use block_array().

Return type:

sparse matrix or array

See also

block_array

Examples

>>> from scipy.sparse import coo_array, bmat
>>> A = coo_array([[1, 2], [3, 4]])
>>> B = coo_array([[5], [6]])
>>> C = coo_array([[7]])
>>> bmat([[A, B], [None, C]]).toarray()
array([[1, 2, 5],
       [3, 4, 6],
       [0, 0, 7]])

>>> bmat([[A, None], [None, C]]).toarray()
array([[1, 2, 0],
       [3, 4, 0],
       [0, 0, 7]])

class common.csr_array(arg1, shape=None, dtype=None, copy=False, *, maxprint=None)

Bases: _csr_base, sparray

Compressed Sparse Row array.

This can be instantiated in several ways:

csr_array(D)

where D is a 2-D ndarray

csr_array(S)

with another sparse array or matrix S (equivalent to S.tocsr())

csr_array((M, N), [dtype])

to construct an empty array with shape (M, N) dtype is optional, defaulting to dtype=’d’.

csr_array((data, (row_ind, col_ind)), [shape=(M, N)])

where data, row_ind and col_ind satisfy the relationship a[row_ind[k], col_ind[k]] = data[k].

csr_array((data, indices, indptr), [shape=(M, N)])

is the standard CSR representation where the column indices for row i are stored in indices[indptr[i]:indptr[i+1]] and their corresponding values are stored in data[indptr[i]:indptr[i+1]]. If the shape parameter is not supplied, the array dimensions are inferred from the index arrays.

dtype

Data type of the array

Type:

dtype

shape

Shape of the array

Type:

2-tuple

ndim

Number of dimensions (this is always 2)

Type:

int

nnz

size

data

CSR format data array of the array

indices

CSR format index array of the array

indptr

CSR format index pointer array of the array

has_sorted_indices

has_canonical_format

T

Notes

Sparse arrays can be used in arithmetic operations: they support addition, subtraction, multiplication, division, and matrix power.

Advantages of the CSR format

• efficient arithmetic operations CSR + CSR, CSR * CSR, etc.

• efficient row slicing

• fast matrix vector products

Disadvantages of the CSR format

• slow column slicing operations (consider CSC)

• changes to the sparsity structure are expensive (consider LIL or DOK)

Canonical Format

• Within each row, indices are sorted by column.

• There are no duplicate entries.

Examples

>>> import numpy as np
>>> from scipy.sparse import csr_array
>>> csr_array((3, 4), dtype=np.int8).toarray()
array([[0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0]], dtype=int8)

>>> row = np.array([0, 0, 1, 2, 2, 2])
>>> col = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> csr_array((data, (row, col)), shape=(3, 3)).toarray()
array([[1, 0, 2],
       [0, 0, 3],
       [4, 5, 6]])

>>> indptr = np.array([0, 2, 3, 6])
>>> indices = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> csr_array((data, indices, indptr), shape=(3, 3)).toarray()
array([[1, 0, 2],
       [0, 0, 3],
       [4, 5, 6]])

Duplicate entries are summed together:

>>> row = np.array([0, 1, 2, 0])
>>> col = np.array([0, 1, 1, 0])
>>> data = np.array([1, 2, 4, 8])
>>> csr_array((data, (row, col)), shape=(3, 3)).toarray()
array([[9, 0, 0],
       [0, 2, 0],
       [0, 4, 0]])

As an example of how to construct a CSR array incrementally, the following snippet builds a term-document array from texts:

>>> docs = [["hello", "world", "hello"], ["goodbye", "cruel", "world"]]
>>> indptr = [0]
>>> indices = []
>>> data = []
>>> vocabulary = {}
>>> for d in docs:
...     for term in d:
...         index = vocabulary.setdefault(term, len(vocabulary))
...         indices.append(index)
...         data.append(1)
...     indptr.append(len(indices))
...
>>> csr_array((data, indices, indptr), dtype=int).toarray()
array([[2, 1, 0, 0],
       [0, 1, 1, 1]])

class common.dia_array(arg1, shape=None, dtype=None, copy=False, *, maxprint=None)

Bases: _dia_base, sparray

Sparse array with DIAgonal storage.

This can be instantiated in several ways:

dia_array(D)

where D is a 2-D ndarray

dia_array(S)

with another sparse array or matrix S (equivalent to S.todia())

dia_array((M, N), [dtype])

to construct an empty array with shape (M, N), dtype is optional, defaulting to dtype=’d’.

dia_array((data, offsets), shape=(M, N))

where the data[k,:] stores the diagonal entries for diagonal offsets[k] (See example below)

dtype

Data type of the array

Type:

dtype

shape

Shape of the array

Type:

2-tuple

ndim

Number of dimensions (this is always 2)

Type:

int

nnz

size

data

DIA format data array of the array

offsets

DIA format offset array of the array

T

Notes

Sparse arrays can be used in arithmetic operations: they support addition, subtraction, multiplication, division, and matrix power. Sparse arrays with DIAgonal storage do not support slicing.

Examples

>>> import numpy as np
>>> from scipy.sparse import dia_array
>>> dia_array((3, 4), dtype=np.int8).toarray()
array([[0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0]], dtype=int8)

>>> data = np.array([[1, 2, 3, 4]]).repeat(3, axis=0)
>>> offsets = np.array([0, -1, 2])
>>> dia_array((data, offsets), shape=(4, 4)).toarray()
array([[1, 0, 3, 0],
       [1, 2, 0, 4],
       [0, 2, 3, 0],
       [0, 0, 3, 4]])

>>> from scipy.sparse import dia_array
>>> n = 10
>>> ex = np.ones(n)
>>> data = np.array([ex, 2 * ex, ex])
>>> offsets = np.array([-1, 0, 1])
>>> dia_array((data, offsets), shape=(n, n)).toarray()
array([[2., 1., 0., ..., 0., 0., 0.],
       [1., 2., 1., ..., 0., 0., 0.],
       [0., 1., 2., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 2., 1., 0.],
       [0., 0., 0., ..., 1., 2., 1.],
       [0., 0., 0., ..., 0., 1., 2.]])

common.spsolve(A, b, permc_spec=None, use_umfpack=True)

Solve the sparse linear system Ax=b, where b may be a vector or a matrix.

Parameters:

• A (ndarray or sparse array or matrix) – The square matrix A will be converted into CSC or CSR form

• b (ndarray or sparse array or matrix) – The matrix or vector representing the right hand side of the equation. If a vector, b.shape must be (n,) or (n, 1).

• permc_spec (str, optional) –

  How to permute the columns of the matrix for sparsity preservation. (default: ‘COLAMD’)

  • NATURAL: natural ordering.

  • MMD_ATA: minimum degree ordering on the structure of A^T A.

  • MMD_AT_PLUS_A: minimum degree ordering on the structure of A^T+A.

  • COLAMD: approximate minimum degree column ordering [1], [2].

• use_umfpack (bool, optional) – if True (default) then use UMFPACK for the solution [3], [4], [5], [6] . This is only referenced if b is a vector and scikits.umfpack is installed.

Returns:

x – the solution of the sparse linear equation. If b is a vector, then x is a vector of size A.shape[1] If b is a matrix, then x is a matrix of size (A.shape[1], b.shape[1])

Return type:

ndarray or sparse array or matrix

Notes

For solving the matrix expression AX = B, this solver assumes the resulting matrix X is sparse, as is often the case for very sparse inputs. If the resulting X is dense, the construction of this sparse result will be relatively expensive. In that case, consider converting A to a dense matrix and using scipy.linalg.solve or its variants.

References

[1]

T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD, an approximate column minimum degree ordering algorithm, ACM Trans. on Mathematical Software, 30(3), 2004, pp. 377–380. :doi:`10.1145/1024074.1024080`

[2]

T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, A column approximate minimum degree ordering algorithm, ACM Trans. on Mathematical Software, 30(3), 2004, pp. 353–376. :doi:`10.1145/1024074.1024079`

[3]

T. A. Davis, Algorithm 832: UMFPACK - an unsymmetric-pattern multifrontal method with a column pre-ordering strategy, ACM Trans. on Mathematical Software, 30(2), 2004, pp. 196–199. https://dl.acm.org/doi/abs/10.1145/992200.992206

[4]

T. A. Davis, A column pre-ordering strategy for the unsymmetric-pattern multifrontal method, ACM Trans. on Mathematical Software, 30(2), 2004, pp. 165–195. https://dl.acm.org/doi/abs/10.1145/992200.992205

[5]

T. A. Davis and I. S. Duff, A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Trans. on Mathematical Software, 25(1), 1999, pp. 1–19. https://doi.org/10.1145/305658.287640

[6]

T. A. Davis and I. S. Duff, An unsymmetric-pattern multifrontal method for sparse LU factorization, SIAM J. Matrix Analysis and Computations, 18(1), 1997, pp. 140–158. https://doi.org/10.1137/S0895479894246905T.

Examples

>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import spsolve
>>> A = csc_array([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
>>> B = csc_array([[2, 0], [-1, 0], [2, 0]], dtype=float)
>>> x = spsolve(A, B)
>>> np.allclose(A.dot(x).toarray(), B.toarray())
True

common.time()

common.vstack(blocks, format=None, dtype=None)

Stack sparse arrays vertically (row wise)

Parameters:

• blocks – sequence of sparse arrays with compatible shapes

• format (str, optional) – sparse format of the result (e.g., “csr”) by default an appropriate sparse array format is returned. This choice is subject to change.

• dtype (dtype, optional) – The data-type of the output array. If not given, the dtype is determined from that of blocks.

Returns:

new_array – If any block in blocks is a sparse array, return a sparse array. Otherwise return a sparse matrix.

If you want a sparse array built from blocks that are not sparse arrays, use block(vstack(blocks)) or convert one block e.g. blocks[0] = csr_array(blocks[0]).

Return type:

sparse matrix or array

See also

hstack

stack sparse matrices horizontally (column wise)

Examples

>>> from scipy.sparse import coo_array, vstack
>>> A = coo_array([[1, 2], [3, 4]])
>>> B = coo_array([[5, 6]])
>>> vstack([A, B]).toarray()
array([[1, 2],
       [3, 4],
       [5, 6]])

config

Configuration file for fuel cell 3D simulation. Contains all configuration parameters including domain dimensions, physical constants, material properties, and simulation settings.

read_image

Image reading module for fuel cell microstructure.

This module provides functionality to read voxelized microstructure data from TIFF image files for solid oxide fuel cell simulations.

read_image.read_in_image(file_name, studied_physics, dimension)

Read a voxelized microstructure image from a TIFF file.

This function loads a 3D TIFF image containing phase labels for the fuel cell microstructure. Each voxel is labeled with a phase ID:

• 0: Pore (gas phase)

• 1: Electrolyte (ionic conductor)

• 2: Electrode (electronic conductor)

Parameters:

• file_name (str) – Path to the TIFF image file.

• studied_physics (str) – Type of physics being studied (“fuel cell” or “battery”).

• dimension (int) – Problem dimension (2 or 3).

Returns:

• img_ (numpy.ndarray) – 3D array of shape (nx, ny, nz) containing phase IDs at each voxel.

• unic_grain_id (list) – List of unique phase IDs found in the image.

• num_pixels_xyz (list) – List [nx, ny, nz] of voxel counts in each direction.

Examples

>>> img, phases, dims = read_in_image("M_3d_3phases_2K.tif", "fuel cell", 3)
>>> print(f"Image shape: {img.shape}")
>>> print(f"Phases found: {phases}")

shape_function_in_domain

RKPM shape function computation module.

This module computes Reproducing Kernel Particle Method (RKPM) shape functions and their derivatives. The RKPM shape function for node I at point x is:

Ψ_I(x) = H^T(0) M^{-1}(x) H(x - x_I) φ(z)

where:

• H(x) is the polynomial basis vector [1, x/h, y/h, z/h]^T

• M(x) is the moment matrix (ensures reproduction of polynomials)

• φ(z) is the cubic B-spline kernel function

• z = ||x - x_I|| / a_I is the normalized distance

The module provides both original loop-based and optimized vectorized implementations for performance.

shape_function_in_domain.compute_phi_M(x_G, Gauss_grain_id, x_nodes, nodes_grain_id, a, M, M_P_x, M_P_y, num_interface_segments, interface_nodes, BxByCxCy, IM_RKPM, single_grain, M_P_z=None, use_vectorized=True, dtype=<class 'numpy.float64'>)

Compute RKPM kernel values φ and moment matrices M at all Gauss points.

This function evaluates the cubic B-spline kernel for all Gauss point-node pairs within the support radius (z ≤ 1) and accumulates the moment matrices needed for shape function computation.

The moment matrix at Gauss point x is:

M(x) = Σ_I H(x - x_I) H^T(x - x_I) φ(z_I)

where the sum is over all nodes I with z_I ≤ 1.

Parameters:

• x_G (np.ndarray) – Gauss point coordinates, shape (n_gauss, dim).

• Gauss_grain_id (np.ndarray) – Phase ID at each Gauss point, shape (n_gauss,).

• x_nodes (np.ndarray) – Node coordinates, shape (n_nodes, dim).

• nodes_grain_id (np.ndarray) – Phase ID at each node, shape (n_nodes,).

• a (np.ndarray) – Support radius for each node, shape (n_nodes,).

• M (np.ndarray) – Moment matrix array to accumulate, shape (n_gauss, basis, basis). For 3D: basis=4 (linear), For 2D: basis=3.

• M_P_x (np.ndarray) – Partial derivatives of M w.r.t. x, y. Same shape as M.

• M_P_y (np.ndarray) – Partial derivatives of M w.r.t. x, y. Same shape as M.

• num_interface_segments (int) – Number of interface segments (for IM-RKPM method).

• interface_nodes (np.ndarray) – Coordinates of nodes on interfaces.

• BxByCxCy (np.ndarray) – Interface segment endpoint coordinates.

• IM_RKPM (str) – “True” for interface-modified RKPM, “False” otherwise.

• single_grain (str) – “True” for single grain, “False” for multi-phase.

• M_P_z (np.ndarray, optional) – Partial derivative of M w.r.t. z (3D only).

• use_vectorized (bool, optional) – If True, use optimized vectorized implementation. Default True.

• dtype (numpy dtype, optional) – Data type for computations. Default np.float64.

Returns:

• phi_nonzero_index_row (np.ndarray) – Gauss point indices for non-zero kernel values, shape (n_nnz,).

• phi_nonzero_index_column (np.ndarray) – Node indices for non-zero kernel values, shape (n_nnz,).

• phi_nonzerovalue_data (np.ndarray) – Kernel values φ(z) at non-zero pairs, shape (n_nnz,).

• phi_P_x_nonzerovalue_data (np.ndarray) – ∂φ/∂x at non-zero pairs, shape (n_nnz,).

• phi_P_y_nonzerovalue_data (np.ndarray) – ∂φ/∂y at non-zero pairs, shape (n_nnz,).

• phi_P_z_nonzerovalue_data (np.ndarray) – ∂φ/∂z at non-zero pairs (3D), shape (n_nnz,).

• M (np.ndarray) – Updated moment matrix, shape (n_gauss, basis, basis).

• M_P_x, M_P_y, M_P_z (np.ndarray) – Updated moment matrix derivatives.

shape_function_in_domain.compute_phi_M_standard(x_G, Gauss_grain_id, x_nodes, nodes_grain_id, a, M, M_P_x, M_P_y, num_interface_segments, interface_nodes, BxByCxCy, M_P_z=None)

Compute phi_M using standard method (else case).

shape_function_in_domain.compute_phi_M_with_interface_method(x_G, Gauss_grain_id, x_nodes, nodes_grain_id, a, M, M_P_x, M_P_y, num_interface_segments, interface_nodes, BxByCxCy, M_P_z=None)

Compute phi_M using interface method (IM_RKPM=True and single_grain=False).

shape_function_in_domain.compute_z_and_H_2d(x_G_i, x_nodes_j, a_j, H_scaling_factor, eps)

Compute z values and H matrices for 2D case (M shape is 3x3).

shape_function_in_domain.compute_z_and_H_3d(x_G_i, x_nodes_j, a_j, H_scaling_factor, eps)

Compute z values and H matrices for 3D case (M shape is 4x4).

shape_function_in_domain.shape_func_n_nodes_by_n_nodes(x_G, x_nodes, num_non_zero_phi_a, HT0, M, phi_nonzerovalue_data, phi_nonzero_index_row, phi_nonzero_index_column, use_vectorized=True)

Compute shape functions for non-zero phi values.

Parameters:

• x_G – Gauss points

• x_nodes – Node points

• num_non_zero_phi_a – Number of non-zero phi values

• HT0 – Initial H transpose matrix

• M – Moment matrices

• phi_nonzerovalue_data – Non-zero phi values

• phi_nonzero_index_row – Row indices for non-zero phi

• phi_nonzero_index_column – Column indices for non-zero phi

• use_vectorized – Whether to use vectorized implementation if available

Returns:

List of shape function values

shape_function_in_domain.shape_func_n_nodes_by_n_nodes_vectorized(x_G, x_nodes, num_non_zero_phi_a, HT0, M, phi_nonzerovalue_data, phi_nonzero_index_row, phi_nonzero_index_column)

Vectorized computation of shape functions for non-zero phi values.

Parameters:

• x_G – Gauss points array (n_gauss, dim)

• x_nodes – Node points array (n_nodes, dim)

• num_non_zero_phi_a – Number of non-zero phi values

• HT0 – Initial H transpose matrix

• M – Moment matrices (n_gauss, basis_size, basis_size)

• phi_nonzerovalue_data – Non-zero phi values

• phi_nonzero_index_row – Row indices for non-zero phi

• phi_nonzero_index_column – Column indices for non-zero phi

Returns:

Array of shape function values of shape (num_non_zero_phi_a,)

Return type:

numpy.ndarray

shape_function_in_domain.shape_grad_shape_func(x_G, x_nodes, num_non_zero_phi_a, HT0, M, M_P_x, M_P_y, differential_method, HT1, HT2, phi_nonzerovalue_data, phi_P_x_nonzerovalue_data, phi_P_y_nonzerovalue_data, phi_nonzero_index_row, phi_nonzero_index_column, det_J_time_weight, IM_RKPM, M_P_z=None, HT3=None, phi_P_z_nonzerovalue_data=None, use_vectorized=True, dtype=<class 'numpy.float64'>)

Compute RKPM shape functions and their gradients.

Evaluates the shape function Ψ_I(x) and its spatial derivatives for all non-zero Gauss-node pairs. The shape function is:

Ψ_I(x) = H^T(0) M^{-1}(x) H(x - x_I) φ(z)

Gradients are computed using the direct differentiation method:

∂Ψ/∂x = H^T M^{-1} H ∂φ/∂x + H^T ∂M^{-1}/∂x H φ + H^T M^{-1} ∂H/∂x φ

Parameters:

• x_G (np.ndarray) – Gauss point coordinates, shape (n_gauss, dim).

• x_nodes (np.ndarray) – Node coordinates, shape (n_nodes, dim).

• num_non_zero_phi_a (int) – Number of non-zero kernel values.

• HT0 (np.ndarray) – Polynomial basis at origin [1, 0, 0, 0], shape (basis,).

• M (np.ndarray) – Moment matrices at Gauss points, shape (n_gauss, basis, basis).

• M_P_x (np.ndarray) – Moment matrix derivatives, same shape as M.

• M_P_y (np.ndarray) – Moment matrix derivatives, same shape as M.

• differential_method (str) – “direct” for direct differentiation, “implicite” for implicit.

• HT1 (np.ndarray) – Basis vectors for implicit method derivatives.

• HT2 (np.ndarray) – Basis vectors for implicit method derivatives.

• phi_nonzerovalue_data (np.ndarray) – Non-zero kernel values, shape (n_nnz,).

• phi_P_*_nonzerovalue_data (np.ndarray) – Kernel derivatives, shape (n_nnz,).

• phi_nonzero_index_row (np.ndarray) – Sparse indices for non-zero pairs, shape (n_nnz,).

• phi_nonzero_index_column (np.ndarray) – Sparse indices for non-zero pairs, shape (n_nnz,).

• det_J_time_weight (np.ndarray) – Jacobian × Gauss weight at each Gauss point, shape (n_gauss,).

• IM_RKPM (str) – “True” for interface-modified RKPM.

• M_P_z (np.ndarray, optional) – Moment matrix z-derivative (3D).

• HT3 (np.ndarray, optional) – Basis vector for implicit z-derivative.

• phi_P_z_nonzerovalue_data (np.ndarray, optional) – Kernel z-derivative values.

• use_vectorized (bool, optional) – Use vectorized implementation. Default True.

• dtype (numpy dtype, optional) – Data type for computations. Default np.float64.

Returns:

• shape_func_value (np.ndarray) – Shape function values, shape (n_nnz,).

• shape_func_times_det_J_time_weight_value (np.ndarray) – Weighted shape functions, shape (n_nnz,).

• grad_shape_func_x_value (np.ndarray) – ∂Ψ/∂x values, shape (n_nnz,).

• grad_shape_func_y_value (np.ndarray) – ∂Ψ/∂y values, shape (n_nnz,).

• grad_shape_func_z_value (np.ndarray) – ∂Ψ/∂z values (3D), shape (n_nnz,).

• grad_shape_func_*_times_det_J_time_weight_value (np.ndarray) – Weighted gradient values, shape (n_nnz,).

shape_function_vectorized

Vectorized version of compute_phi_M_standard and related functions.

This module provides vectorized implementations that operate on all Gauss points and nodes simultaneously, eliminating the need for explicit loops.

shape_function_vectorized.compute_phi_M_standard_sparse(x_G, Gauss_grain_id, x_nodes, nodes_grain_id, a, M, M_P_x, M_P_y, num_interface_segments, interface_nodes, BxByCxCy, M_P_z=None, dtype=<class 'numpy.float64'>)

Optimized sparse version of compute_phi_M_standard using vectorized operations.

This function eliminates the loop over Gauss points by using advanced indexing and np.add.at for efficient accumulation. It provides 10-100x speedup for large numbers of Gauss points.

Parameters:

dtype – Data type for computations (default: np.float64).

shape_function_vectorized.compute_phi_M_standard_vectorized(x_G, Gauss_grain_id, x_nodes, nodes_grain_id, a, M, M_P_x, M_P_y, num_interface_segments, interface_nodes, BxByCxCy, M_P_z=None, dtype=<class 'numpy.float64'>)

Vectorized version of compute_phi_M_standard.

This function computes all Gauss-node interactions simultaneously using NumPy broadcasting and vectorized operations.

Parameters:

dtype – Data type for computations (default: np.float32 for memory efficiency). Use np.float64 for higher precision if needed.

shape_function_vectorized.compute_phi_kernel(z)

Compute phi kernel values based on distance z.

The kernel is a cubic spline-like function: - For 0 <= z < 0.5: phi = 2/3 - 4*z^2 + 4*z^3 - For 0.5 <= z <= 1: phi = 4/3 - 4*z + 4*z^2 - 4/3*z^3 - For z > 1: phi = 0

Parameters:

z – Array of normalized distances

Returns:

Array of kernel values phi_P_z: Array of kernel derivatives w.r.t. z

Return type:

phi

shape_function_vectorized.compute_z_and_H_2d_vectorized(x_G, x_nodes, a, H_scaling_factor, eps, dtype=<class 'numpy.float64'>)

Compute z values and H matrices for 2D case (vectorized).

Parameters:

• x_G – Array of Gauss points, shape (n_gauss, 2)

• x_nodes – Array of node points, shape (n_nodes, 2)

• a – Array of scaling factors, shape (n_nodes,)

• H_scaling_factor – Scaling factor for H matrices

• eps – Small value to avoid division by zero

• dtype – Data type for arrays (default: np.float32 for memory efficiency)

Returns:

Tuple containing vectorized z values, derivatives, and H matrices

shape_function_vectorized.compute_z_and_H_3d_vectorized(x_G, x_nodes, a, H_scaling_factor, eps, dtype=<class 'numpy.float64'>)

Compute z values and H matrices for 3D case (vectorized).

Parameters:

• x_G – Array of Gauss points, shape (n_gauss, 3)

• x_nodes – Array of node points, shape (n_nodes, 3)

• a – Array of scaling factors, shape (n_nodes,)

• H_scaling_factor – Scaling factor for H matrices

• eps – Small value to avoid division by zero

• dtype – Data type for arrays (default: np.float32 for memory efficiency)

Returns:

Tuple containing vectorized z values, derivatives, and H matrices

shape_grad_func_vectorized

Vectorized version of shape_grad_shape_func and related functions.

This module provides vectorized implementations for computing shape functions and their gradients, with optimized batch matrix operations.

shape_grad_func_vectorized.batch_matrix_inverse(matrices)

Compute inverses of multiple matrices in a batch.

Parameters:

matrices – Array of shape (n_matrices, dim, dim)

Returns:

Array of inverted matrices with same shape

shape_grad_func_vectorized.compute_H_matrices_sparse_2d(x_G, x_nodes, gauss_indices, node_indices, H_scaling_factor, dtype=<class 'numpy.float64'>)

Compute H matrices for sparse Gauss-node pairs in 2D.

Parameters:

• x_G – Gauss points array

• x_nodes – Node points array

• gauss_indices – Indices of Gauss points in sparse pairs

• node_indices – Indices of nodes in sparse pairs

• H_scaling_factor – Scaling factor for H matrices

• dtype – Data type for computations

Returns:

H_T, HT_P_x, HT_P_y, H, H_P_x, H_P_y for all sparse pairs

shape_grad_func_vectorized.compute_H_matrices_sparse_3d(x_G, x_nodes, gauss_indices, node_indices, H_scaling_factor, dtype=<class 'numpy.float64'>)

Compute H matrices for sparse Gauss-node pairs in 3D.

Parameters:

• x_G – Gauss points array

• x_nodes – Node points array

• gauss_indices – Indices of Gauss points in sparse pairs

• node_indices – Indices of nodes in sparse pairs

• H_scaling_factor – Scaling factor for H matrices

• dtype – Data type for computations

Returns:

H_T, HT_P_x, HT_P_y, HT_P_z, H, H_P_x, H_P_y, H_P_z for all sparse pairs

shape_grad_func_vectorized.shape_grad_shape_func_vectorized(x_G, x_nodes, num_non_zero_phi_a, HT0, M, M_P_x, M_P_y, differential_method, HT1, HT2, phi_nonzerovalue_data, phi_P_x_nonzerovalue_data, phi_P_y_nonzerovalue_data, phi_nonzero_index_row, phi_nonzero_index_column, det_J_time_weight, IM_RKPM, M_P_z=None, HT3=None, phi_P_z_nonzerovalue_data=None, dtype=<class 'numpy.float64'>)

Vectorized version of shape_grad_shape_func.

This function computes shape functions and gradients using vectorized operations and batch matrix inversions for improved performance.

Parameters:

dtype – Data type for computations (np.float32 or np.float64)

get_nodes_gauss_points_vectorized

Fully vectorized versions of Gauss point computation functions. This module provides high-performance implementations using NumPy broadcasting and vectorization to eliminate loops.

get_nodes_gauss_points_vectorized.x_G_and_def_J_time_weight_3d_fuelcell_domain_vectorized(cell_nodes_x, cell_nodes_y, cell_nodes_z, x_G_domain, weight_G_domain)

Fully vectorized computation of Gauss points and Jacobian determinants for 3D hexahedral cells.

Parameters:

cell_nodes_xnp.ndarray

Array of shape (n_cells, 8) containing x-coordinates of cell vertices

cell_nodes_ynp.ndarray

Array of shape (n_cells, 8) containing y-coordinates of cell vertices

cell_nodes_znp.ndarray

Array of shape (n_cells, 8) containing z-coordinates of cell vertices

x_G_domainnp.ndarray or list

Array of shape (n_gauss, 3) containing Gauss points in reference coordinates

weight_G_domainnp.ndarray or list

Array of shape (n_gauss,) containing Gauss weights

Returns:

x_Gnp.ndarray

Array of shape (n_cells * n_gauss, 3) containing [x, y, z] coordinates of Gauss points

det_J_time_weightnp.ndarray

Array of shape (n_cells * n_gauss,) containing Jacobian determinants times weights

get_nodes_gauss_points_vectorized.x_G_b_and_det_J_b_time_weight_3d_fuelcell_2d_boundary_interface_vectorized(cell_nodes_boundary_x, cell_nodes_boundary_y, cell_nodes_boundary_z, x_G_domain, weight_G_domain)

Fully vectorized computation of Gauss points and Jacobian determinants for 2D boundary interfaces in a 3D fuel cell domain. This version processes all cells of each type simultaneously.

Parameters:

cell_nodes_boundary_x/y/znp.ndarray

Arrays of shape (n_cells, 4) containing coordinates of vertices for each boundary cell

x_G_domainnp.ndarray

Array of shape (n_gauss, 2) containing reference Gauss quadrature points

weight_G_domainnp.ndarray

Array of shape (n_gauss,) containing Gauss weights

Returns:

x_Glist

List of [x, y, z] coordinates for all Gauss points

det_J_time_weightlist

List of Jacobian determinants multiplied by weights

get_nodes_gauss_points_vectorized.x_G_b_and_det_J_b_time_weight_3d_fuelcell_2d_boundary_vectorized(cell_nodes_boundary_x, cell_nodes_boundary_z, y_coords_on_boundary, x_G_domain, weight_G_domain)

Fully vectorized computation of Gauss points and Jacobian determinants for 2D boundaries in a 3D fuel cell domain with fixed y-coordinate.

Parameters:

cell_nodes_boundary_xnp.ndarray

Array of shape (n_cells, 4) containing x-coordinates of boundary cell vertices

cell_nodes_boundary_znp.ndarray

Array of shape (n_cells, 4) containing z-coordinates of boundary cell vertices

y_coords_on_boundaryfloat

Fixed y-coordinate value for all points on this boundary

x_G_domainnp.ndarray

Array of shape (n_gauss, 2) containing reference Gauss quadrature points

weight_G_domainnp.ndarray

Array of shape (n_gauss,) containing Gauss weights

Returns:

x_Glist

List of [x, y, z] coordinates for all Gauss points

det_J_time_weightlist

List of Jacobian determinants multiplied by weights

define_diffusion_matrix_form

Diffusion matrix assembly module using Nitsche’s method.

This module assembles the stiffness matrix and force vector for the steady-state diffusion equation with Dirichlet boundary conditions enforced using Nitsche’s method (weak enforcement).

The discretized system is: K @ c = f

where K = K1 + K2 + K3:

• K1: Domain diffusion integral

• K2: Flux consistency term on boundary

• K3: Nitsche penalty term on boundary

define_diffusion_matrix_form.diffusion_matrix_fuel_cell(dimension: int, point_or_line_source: ndarray, shape_func_point_or_line_nodes: csr_array, g_diretchlet: ndarray, beta_Nitsche: ndarray, normal_vector_x: ndarray, normal_vector_y: ndarray, global_diffusion: ndarray, grad_shape_func_x: csr_array, grad_shape_func_y: csr_array, grad_shape_func_x_times_det_J_time_weight: csr_array, grad_shape_func_y_times_det_J_time_weight: csr_array, shape_func_b: csr_array, shape_func_b_times_det_J_b_time_weight: csr_array, grad_shape_func_b_x_times_det_J_b_time_weight: csr_array, grad_shape_func_b_y_times_det_J_b_time_weight: csr_array, shape_func_inter_times_det_J_b_time_weight: csr_array | None = None, interface_source: ndarray | None = None, grad_shape_func_z: csr_array | None = None, grad_shape_func_z_times_det_J_time_weight: csr_array | None = None, grad_shape_func_b_z_times_det_J_b_time_weight: csr_array | None = None, normal_vector_z: ndarray | None = None) → Tuple[ndarray, ndarray]

Assemble diffusion stiffness matrix and force vector using Nitsche’s method.

This function constructs the linear system for the steady-state diffusion equation with Dirichlet boundary conditions enforced weakly via Nitsche’s method. Used for delta function (point) source terms.

Parameters:

• dimension (int) – Problem dimension (2 or 3).

• point_or_line_source (np.ndarray) – Source term values at source points, shape (n_source,).

• shape_func_point_or_line_nodes (csr_array) – Shape functions at source points, shape (n_source, n_nodes).

• g_diretchlet (np.ndarray) – Dirichlet boundary values, shape (n_gauss_boundary,).

• beta_Nitsche (np.ndarray) – Nitsche penalty parameter, shape (n_gauss_boundary,).

• normal_vector_x (np.ndarray) – Outward normal components at boundary Gauss points.

• normal_vector_y (np.ndarray) – Outward normal components at boundary Gauss points.

• global_diffusion (np.ndarray) – Diffusion coefficient at domain Gauss points.

• grad_shape_func_x (csr_array) – Shape function gradients, shape (n_gauss, n_nodes).

• grad_shape_func_y (csr_array) – Shape function gradients, shape (n_gauss, n_nodes).

• grad_shape_func_x_times_det_J_time_weight (csr_array) – Weighted gradient matrices for integration, shape (n_gauss, n_nodes).

• shape_func_b (csr_array) – Shape functions at boundary Gauss points.

• shape_func_b_times_det_J_b_time_weight (csr_array) – Weighted boundary shape functions.

• grad_shape_func_b_*_times_det_J_b_time_weight (csr_array) – Weighted boundary gradient matrices.

• shape_func_inter_times_det_J_b_time_weight (csr_array, optional) – Shape functions at interface for interface source.

• interface_source (np.ndarray, optional) – Interface source term values.

• grad_shape_func_z (csr_array, optional) – z-gradient of shape functions (3D only).

• normal_vector_z (np.ndarray, optional) – z-component of outward normal (3D only).

Returns:

• K (np.ndarray) – Stiffness matrix, shape (n_nodes, n_nodes).

• f (np.ndarray) – Force vector, shape (n_nodes,).

define_diffusion_matrix_form.diffusion_matrix_fuel_cell_distributed_point_source(dimension: int, distributed_point_or_line_source: ndarray, shape_func_distributed_point_or_line_nodes: csr_array, g_diretchlet: ndarray, beta_Nitsche: ndarray, normal_vector_x: ndarray, normal_vector_y: ndarray, global_diffusion: ndarray, grad_shape_func_x: csr_array, grad_shape_func_y: csr_array, grad_shape_func_x_times_det_J_time_weight: csr_array, grad_shape_func_y_times_det_J_time_weight: csr_array, shape_func_b: csr_array, shape_func_b_times_det_J_b_time_weight: csr_array, grad_shape_func_b_x_times_det_J_b_time_weight: csr_array, grad_shape_func_b_y_times_det_J_b_time_weight: csr_array, shape_func_inter_times_det_J_b_time_weight: csr_array | None = None, interface_source: ndarray | None = None, grad_shape_func_z: csr_array | None = None, grad_shape_func_z_times_det_J_time_weight: csr_array | None = None, grad_shape_func_b_z_times_det_J_b_time_weight: csr_array | None = None, normal_vector_z: ndarray | None = None) → Tuple[ndarray, ndarray]

define_mechanical_stiffness_matrix

Mechanical stiffness matrix assembly for 3D linear elasticity.

This module assembles the mechanical stiffness matrix for 3D elasticity problems in solid oxide fuel cells. It includes:

• Construction of the 6x6 elasticity tensor (Voigt notation)

• Rotation for anisotropic materials

• Damage model degradation

• Full 3x3 block stiffness matrix assembly

The mechanical stiffness matrix relates nodal displacements to forces:

K_mech @ u = f

where u = [u_x, u_y, u_z] and K_mech is a 3x3 block matrix.

define_mechanical_stiffness_matrix.mechanical_C_tensor_3d(num_gauss_points_in_domain: int, D_damage: ndarray, lambda_mechanical: ndarray, mu: ndarray, gauss_angle: ndarray, gauss_rotation_axis: ndarray) → Tuple[ndarray, ndarray]

Construct the 3D elasticity tensor with rotation and damage.

Computes the 6x6 elasticity tensor in Voigt notation for each Gauss point, including material rotation and damage degradation.

The elasticity tensor relates stress to strain:

σ = C : (ε - ε*)

where ε* is the chemical expansion strain.

Parameters:

• num_gauss_points_in_domain (int) – Number of Gauss points in the domain.

• D_damage (np.ndarray) – Damage variable at each Gauss point, shape (n_gauss, 1). Value 0 = undamaged, 1 = fully damaged.

• lambda_mechanical (np.ndarray) – First Lamé constant (λ) at each Gauss point.

• mu (np.ndarray) – Second Lamé constant (shear modulus, μ) at each Gauss point.

• gauss_angle (np.ndarray) – Rotation angle at each Gauss point for anisotropic materials.

• gauss_rotation_axis (np.ndarray) – Rotation axis vectors, shape (n_gauss, 3).

Returns:

• C (np.ndarray) – Elasticity tensor, shape (6, 6, n_gauss, 1). Components are in Voigt notation: [11, 22, 33, 23, 13, 12].

• T_c (np.ndarray) – Rotation transformation matrix, shape (6, 6, n_gauss).

define_mechanical_stiffness_matrix.mechanical_force_matrix_3d(x_G, C, epsilon_D1, epsilon_D2, epsilon_D3, epsilon_D4, epsilon_D5, epsilon_D6, grad_shape_func_x_times_det_J_time_weight, grad_shape_func_y_times_det_J_time_weight, grad_shape_func_z_times_det_J_time_weight)

define_mechanical_stiffness_matrix.mechanical_stiffness_matrix_3d_fuel_cell(C: ndarray, num_gauss_points_in_domain: int, grad_shape_func_x_times_det_J_time_weight: csr_array, grad_shape_func_x: csr_array, grad_shape_func_y_times_det_J_time_weight: csr_array, grad_shape_func_y: csr_array, grad_shape_func_z_times_det_J_time_weight: csr_array, grad_shape_func_z: csr_array, beta_Nitsche: ndarray, shape_func_fixed_point: csr_array, shape_func_times_det_J_time_weight_fixed_point: csr_array, grad_shape_func_x_fixed_point: csr_array, grad_shape_func_x_times_det_J_time_weight_fixed_point: csr_array, grad_shape_func_y_fixed_point: csr_array, grad_shape_func_y_times_det_J_time_weight_fixed_point: csr_array, grad_shape_func_z_fixed_point: csr_array, grad_shape_func_z_times_det_J_time_weight_fixed_point: csr_array, normal_vector_x_electrolyte: ndarray, normal_vector_y_electrolyte: ndarray, normal_vector_z_electrolyte: ndarray, shape_func_b_electrolyte: csr_array, shape_func_b_times_det_J_b_time_weight_electrolyte: csr_array, grad_shape_func_b_x_electrolyte: csr_array, grad_shape_func_b_x_times_det_J_b_time_weight_electrolyte: csr_array, grad_shape_func_b_y_electrolyte: csr_array, grad_shape_func_b_y_times_det_J_b_time_weight_electrolyte: csr_array, grad_shape_func_b_z_electrolyte: csr_array, grad_shape_func_b_z_times_det_J_b_time_weight_electrolyte: csr_array) → csr_array

Assemble the 3D mechanical stiffness matrix for fuel cell simulation.

Constructs the full 3×3 block stiffness matrix for 3D linear elasticity using RKPM shape functions. Includes both domain integrals and boundary terms (Nitsche’s method for fixed displacement BCs).

The stiffness matrix structure:

K_mech = | K11 K12 K13 |

K21 K22 K23 |

K31 K32 K33 |

where each Kij block is (n_nodes × n_nodes) and couples displacement components i and j.

Parameters:

• C (np.ndarray) – Elasticity tensor, shape (6, 6, n_gauss, 1).

• num_gauss_points_in_domain (int) – Number of Gauss points in the domain.

• grad_shape_func_x (csr_array) – Shape function gradients in domain, each (n_gauss, n_nodes).

• grad_shape_func_y (csr_array) – Shape function gradients in domain, each (n_gauss, n_nodes).

• grad_shape_func_z (csr_array) – Shape function gradients in domain, each (n_gauss, n_nodes).

• grad_shape_func_*_times_det_J_time_weight (csr_array) – Weighted gradient matrices for integration.

• beta_Nitsche (np.ndarray) – Nitsche penalty parameter at boundary points.

• shape_func_fixed_point (csr_array) – Shape functions at fixed displacement points.

• normal_vector_*_electrolyte (np.ndarray) – Normal vector components at boundary Gauss points.

• shape_func_b_electrolyte (csr_array) – Shape functions at boundary Gauss points.

• grad_shape_func_b_*_electrolyte (csr_array) – Gradient matrices at boundary Gauss points.

Returns:

K_mechanical – Mechanical stiffness matrix, shape (3*n_nodes, 3*n_nodes). Ordering: [u_x nodes, u_y nodes, u_z nodes].

Return type:

csr_array