fom.fem_utils
Utilities for finite element (FE) workflows—designed to simplify mesh handling, accelerate solvers, and streamline nonlinear problem setups.
Summary: A collection of helper functions covering mesh and basis manipulation, fast linear and nonlinear solvers, domain loading, and attribute cleanup. These utilities reduce boilerplate and improve performance across FEM pipelines.
Purpose: To provide reusable, composable building blocks that enhance efficiency and code clarity in finite element simulations—particularly when working with iterative solvers, region-partitioned computations, and domain abstractions.
Author: Suparno Bhattacharyya
Functions
| Name | Description |
|---|---|
| build_pc_amgsa | Build a fast algebraic multigrid (AMG) preconditioner. |
| compute_basis_regions | Restrict basis functions to specific mesh regions using element masks. |
| element2location | Retrieve spatial coordinates for each element in a mesh. |
| load_domain | Load all outputs from domain() as attributes on an object. |
| load_mesh_and_basis | Load only mesh and basis from a domain definition. |
| newton_solver | Newton nonlinear solver with optional PETSc acceleration. |
| petsc_solve_csr | Solve A x = b using PETSc KSP with CSR sparse matrix input. |
| unwrap_attr | Unwrap a 0-dimensional NumPy object array attribute to its plain value. |
build_pc_amgsa
fom.fem_utils.build_pc_amgsa(A, **kwargs)Constructs an algebraic multigrid (AMG) preconditioner to accelerate convergence of large-scale iterative solvers, such as Conjugate Gradient (CG). AMG preconditioners are especially effective for elliptic PDE systems arising in structural and thermal FEM problems.
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| A | scipy.sparse matrix or array_like | System matrix to precondition. | required |
| **kwargs | dict | Additional options forwarded to the AMG solver. | {} |
Returns
| Name | Type | Description |
|---|---|---|
| M | scipy.sparse.linalg.LinearOperator | AMG preconditioner compatible with iterative solvers. |
compute_basis_regions
fom.fem_utils.compute_basis_regions(basis, masks)Partitions a global basis into multiple region-specific bases using boolean masks. Each mask identifies which elements belong to a given region, enabling localized assembly or region-wise computations.
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| basis | object | Basis object exposing nelems and a with_elements method. |
required |
| masks | dict of str to ndarray of bool | Maps region names to boolean arrays selecting the corresponding elements. | required |
Returns
| Name | Type | Description |
|---|---|---|
| region_bases | dict of str to object | Maps each region name to its corresponding restricted basis object. |
element2location
fom.fem_utils.element2location(mesh)Computes the spatial coordinates associated with each element by indexing mesh node positions through the connectivity array. Useful for element-level geometric computations such as quadrature point mapping or centroid evaluation.
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| mesh | object | Mesh object with p (node coordinates, shape [spatial_dim, n_nodes]) and t (connectivity, shape [n_nodes_per_element, n_elements]). |
required |
Returns
| Name | Type | Description |
|---|---|---|
| element_coords | ndarray of shape (n_elements, n_local_nodes) | Per-element node coordinates; each row corresponds to one element. |
load_domain
fom.fem_utils.load_domain(instance)Calls the domain() method on the given object and assigns all key-value pairs in the returned dictionary as attributes directly on the object.
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| instance | object | Any object implementing a domain() method that returns a dictionary of attributes. |
required |
Notes
Modifies the instance in-place. After calling this function, all domain fields are accessible directly as instance.attr—eliminating manual unpacking.
load_mesh_and_basis
fom.fem_utils.load_mesh_and_basis(instance)A focused variant of load_domain that loads only the mesh and basis fields from a domain() output, ignoring all other returned data.
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| instance | object | Object whose domain() method returns at least mesh and basis keys. |
required |
Notes
Only mesh and basis are attached to the instance; all other fields are discarded. Instance is modified in-place.
newton_solver
fom.fem_utils.newton_solver(
assemble_fn,
rhs_fn,
u0,
dirichlet_dofs=None,
dirichlet_vals=None,
*assemble_args,
rhs_args=(),
tol=0.01,
maxit=50,
alpha=1.0,
jac_conditioner=False,
ksp_type='cg',
pc_type='ilu',
ksp_rtol=1e-08,
ksp_atol=0.0,
ksp_max_it=2000,
petsc_options=None,
reuse_ksp=True,
force_backend='auto',
)Solves nonlinear systems of equations using Newton’s method. At each iteration, the Jacobian is assembled and the resulting linear system is solved using PETSc KSP (if available) or a NumPy/SciPy fallback. Supports Dirichlet boundary condition enforcement, damped Newton updates via alpha, and optional reuse of the Krylov solver across iterations.
Key Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| assemble_fn | callable | Assembles the Jacobian (tangent stiffness matrix). | required |
| rhs_fn | callable | Assembles the residual (right-hand side) vector. | required |
| u0 | ndarray | Initial guess for the solution. | required |
| dirichlet_dofs | array-like or None | Indices of Dirichlet-constrained degrees of freedom. | None |
| dirichlet_vals | array-like or None | Prescribed values at Dirichlet DOFs. | None |
| tol | float | Convergence tolerance on the residual norm. | 0.01 |
| maxit | int | Maximum number of Newton iterations. | 50 |
| alpha | float | Damping factor for Newton update step. | 1.0 |
| jac_conditioner | bool | Apply conditioning to the Jacobian before solving. | False |
| ksp_type | str | Krylov solver type (e.g., 'cg', 'gmres'). |
'cg' |
| pc_type | str | Preconditioner type (e.g., 'ilu', 'jacobi'). |
'ilu' |
| ksp_rtol | float | Relative tolerance for the Krylov solver. | 1e-08 |
| ksp_atol | float | Absolute tolerance for the Krylov solver. | 0.0 |
| ksp_max_it | int | Maximum Krylov iterations per Newton step. | 2000 |
| petsc_options | dict or None | Additional PETSc runtime options. | None |
| reuse_ksp | bool | Reuse the KSP object across iterations for efficiency. | True |
| force_backend | str | Force solver backend: 'auto', 'petsc', or 'scipy'. |
'auto' |
petsc_solve_csr
fom.fem_utils.petsc_solve_csr(
A_csr,
b,
*,
ksp_type='cg',
pc_type='ilu',
rtol=1e-08,
atol=0.0,
max_it=2000,
petsc_options=None,
cache=None,
)Solves the linear system ( Ax = b ) using PETSc’s Krylov subspace solvers. Accepts a SciPy CSR matrix and converts it internally to PETSc format. An optional solver cache allows reuse of KSP objects, significantly reducing overhead in repeated solves (e.g., inside a Newton loop).
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| A_csr | scipy.sparse.csr_matrix | System matrix in CSR format. | required |
| b | ndarray | Right-hand side vector. | required |
| ksp_type | str | Krylov method (e.g., 'cg', 'gmres', 'bicg'). |
'cg' |
| pc_type | str | Preconditioner (e.g., 'ilu', 'jacobi', 'none'). |
'ilu' |
| rtol | float | Relative convergence tolerance. | 1e-08 |
| atol | float | Absolute convergence tolerance. | 0.0 |
| max_it | int | Maximum number of solver iterations. | 2000 |
| petsc_options | dict or None | Additional PETSc runtime configuration. | None |
| cache | dict or None | Solver cache for KSP reuse across calls. | None |
unwrap_attr
fom.fem_utils.unwrap_attr(instance, attr_name)Inspects a named attribute on an object and, if it is a 0-dimensional NumPy array with object dtype, replaces it with the underlying scalar or Python object it wraps. This commonly arises when loading data from .mat files or NumPy archives where scalars are inadvertently boxed into 0-d arrays.
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| instance | object | Object holding the attribute to unwrap. | required |
| attr_name | str | Name of the attribute to inspect and unwrap. | required |
Notes
Only takes effect when the attribute is a 0-dimensional NumPy array with dtype=object. The instance is modified in-place; after the call, the attribute holds a plain Python value rather than a NumPy wrapper.