CPU Backend

The CPU backend provides all core functionality for topology optimization using CPU-based computations.

CPU backend public API.

Importing from this module gives access to CPU implementations of core pyFANTOM components (meshes, filters, kernels, solvers, optimizers and finite-element helpers). Users typically switch to the CPU backend with:

>>> from pyFANTOM.CPU import StructuredMesh2D, FiniteElement, MinimumCompliance

Mesh Classes

Filter Classes

Finite Element Classes

Stiffness Kernel Classes

Solver Classes

Optimizer Classes

Problem Classes

Detailed Documentation

class pyFANTOM.geom.CPU._mesh.StructuredMesh2D(nx, ny, lx, ly, dtype=<class 'numpy.float64'>, physics: ~pyFANTOM.physics._physx.Physx = <pyFANTOM.physics.LinearElasticity.LinearElasticity object>)[source]

Bases: StructuredMesh

2D structured mesh with uniform rectangular elements.

Creates a regular grid of quadrilateral elements for 2D topology optimization and finite element analysis. All elements have identical geometry, enabling optimized assembly using a single precomputed stiffness matrix.

Parameters:

  • nx (int) – Number of elements in x-direction

  • ny (int) – Number of elements in y-direction

  • lx (float) – Physical length of domain in x-direction

  • ly (float) – Physical length of domain in y-direction

  • dtype (np.dtype, optional) – Data type for arrays (default: np.float64)

  • physics (Physx, optional) – Physics model defining material behavior and element formulation (default: LinearElasticity(E=1.0, nu=1/3))

nelx, nely

Number of elements in x and y directions

Type:

int

dx, dy

Element dimensions

Type:

float

elements

Element connectivity array, shape (nx*ny, 4) with node indices

Type:

ndarray

nodes

Node coordinates, shape (n_nodes, 2)

Type:

ndarray

K_single

Precomputed element stiffness matrix for the uniform geometry

Type:

ndarray

As

Element areas/volumes (all identical for structured mesh)

Type:

ndarray

volume

Total domain volume

Type:

float

dof

Degrees of freedom per node (2 for 2D elasticity)

Type:

int

centroids

Element centroid coordinates, shape (nx*ny, 2)

Type:

ndarray

Notes

  • Element node ordering is counter-clockwise starting from bottom-left

  • Negative areas indicate incorrect node ordering

  • Structured meshes enable fast node-basis assembly via StructuredStiffnessKernel

Examples

>>> from pyFANTOM.CPU import StructuredMesh2D
>>> from pyFANTOM import LinearElasticity
>>> physics = LinearElasticity(E=1.0, nu=0.3, type='PlaneStress')
>>> mesh = StructuredMesh2D(nx=64, ny=32, lx=2.0, ly=1.0, physics=physics)
>>> print(f"Elements: {len(mesh.elements)}, Nodes: {len(mesh.nodes)}")
__init__(nx, ny, lx, ly, dtype=<class 'numpy.float64'>, physics: ~pyFANTOM.physics._physx.Physx = <pyFANTOM.physics.LinearElasticity.LinearElasticity object>)[source]
class pyFANTOM.geom.CPU._mesh.StructuredMesh3D(nx, ny, nz, lx, ly, lz, dtype=<class 'numpy.float64'>, physics: ~pyFANTOM.physics._physx.Physx = <pyFANTOM.physics.LinearElasticity.LinearElasticity object>)[source]

Bases: StructuredMesh

3D structured mesh with uniform hexahedral elements.

Creates a regular grid of hexahedral (brick) elements for 3D topology optimization and finite element analysis. All elements have identical geometry, enabling optimized assembly using a single precomputed stiffness matrix.

Parameters:

  • nx (int) – Number of elements in x-direction

  • ny (int) – Number of elements in y-direction

  • nz (int) – Number of elements in z-direction

  • lx (float) – Physical length of domain in x-direction

  • ly (float) – Physical length of domain in y-direction

  • lz (float) – Physical length of domain in z-direction

  • dtype (np.dtype, optional) – Data type for arrays (default: np.float64)

  • physics (Physx, optional) – Physics model defining material behavior and element formulation (default: LinearElasticity(E=1.0, nu=1/3))

nelx, nely, nelz

Number of elements in x, y, and z directions

Type:

int

dx, dy, dz

Element dimensions

Type:

float

elements

Element connectivity array, shape (nx*ny*nz, 8) with node indices

Type:

ndarray

nodes

Node coordinates, shape (n_nodes, 3)

Type:

ndarray

K_single

Precomputed element stiffness matrix for the uniform geometry

Type:

ndarray

As

Element volumes (all identical for structured mesh)

Type:

ndarray

volume

Total domain volume

Type:

float

dof

Degrees of freedom per node (3 for 3D elasticity)

Type:

int

centroids

Element centroid coordinates, shape (nx*ny*nz, 3)

Type:

ndarray

Notes

  • Element node ordering: bottom face counter-clockwise, then top face counter-clockwise

  • Negative volumes indicate incorrect node ordering

  • Structured 3D meshes are memory-intensive; use MultiGrid solver for large problems

Examples

>>> from pyFANTOM.CPU import StructuredMesh3D
>>> from pyFANTOM import LinearElasticity
>>> physics = LinearElasticity(E=1.0, nu=0.3)
>>> mesh = StructuredMesh3D(nx=32, ny=32, nz=16, lx=1.0, ly=1.0, lz=0.5, physics=physics)
>>> print(f"Elements: {len(mesh.elements)}, DOF: {len(mesh.nodes)*mesh.dof}")
__init__(nx, ny, nz, lx, ly, lz, dtype=<class 'numpy.float64'>, physics: ~pyFANTOM.physics._physx.Physx = <pyFANTOM.physics.LinearElasticity.LinearElasticity object>)[source]
class pyFANTOM.geom.CPU._mesh.GeneralMesh(nodes, elements, dtype=<class 'numpy.float64'>, physics: ~pyFANTOM.physics._physx.Physx = <pyFANTOM.physics.LinearElasticity.LinearElasticity object>)[source]

Bases: object

General unstructured mesh supporting arbitrary element topologies.

Handles unstructured meshes with heterogeneous element types and sizes. Automatically detects whether all elements are uniform (same number of nodes) and stores stiffness matrices efficiently. Removes redundant nodes during initialization.

Parameters:

  • nodes (ndarray) – Node coordinates, shape (n_nodes, spatial_dim)

  • elements (list or ndarray) – Element connectivity. For uniform meshes: array of shape (n_elements, nodes_per_element). For heterogeneous meshes: list of arrays with varying lengths

  • dtype (np.dtype, optional) – Data type for arrays (default: np.float64)

  • physics (Physx, optional) – Physics model defining material behavior and element formulation (default: LinearElasticity(E=1.0, nu=1/3))

nodes

Node coordinates (cleaned of redundant nodes)

Type:

ndarray

elements

Element connectivity arrays

Type:

list or ndarray

is_uniform

True if all elements have the same number of nodes

Type:

bool

elements_flat

Flattened element connectivity for efficient access

Type:

ndarray

elements_ptr

Pointer array for accessing elements in flat storage (for heterogeneous meshes)

Type:

ndarray

Ks

Element stiffness matrices (uniform case)

Type:

ndarray

K_flat

Flattened stiffness matrices (heterogeneous case)

Type:

ndarray

K_ptr

Pointer array for accessing stiffness matrices (heterogeneous case)

Type:

ndarray

As

Element areas/volumes, shape (n_elements,)

Type:

ndarray

volume

Total domain volume

Type:

float

dof

Degrees of freedom per node

Type:

int

centroids

Element centroid coordinates

Type:

ndarray

Notes

  • Automatically removes unused nodes and remaps element connectivity

  • For heterogeneous meshes, uses flat storage with pointer arrays for memory efficiency

  • Use GeneralStiffnessKernel for assembly with general meshes

Examples

>>> from pyFANTOM.CPU import GeneralMesh
>>> from pyFANTOM import LinearElasticity
>>> # Triangular mesh
>>> nodes = np.array([[0,0], [1,0], [0,1], [1,1]])
>>> elements = np.array([[0,1,2], [1,3,2]])
>>> mesh = GeneralMesh(nodes, elements, physics=LinearElasticity())
>>> print(f"Uniform mesh: {mesh.is_uniform}, Elements: {len(mesh.elements)}")
__init__(nodes, elements, dtype=<class 'numpy.float64'>, physics: ~pyFANTOM.physics._physx.Physx = <pyFANTOM.physics.LinearElasticity.LinearElasticity object>)[source]
class pyFANTOM.geom.CPU._filters.StructuredFilter2D(mesh: StructuredMesh2D, r_min)[source]

Bases: FilterKernel

Convolution-based density filter for 2D structured meshes.

Implements smoothing filter for 2D topology optimization using efficient convolution. Ensures minimum feature sizes and mesh-independent solutions.

Parameters:
shape

Filter matrix dimensions (n_elements, n_elements)

Type:

tuple

weights

Precomputed filter weights

Type:

ndarray

offsets

Element index offsets for neighbors

Type:

ndarray

normalizer

Per-element normalization for adjoint

Type:

ndarray

dot(rho)

Forward filter application

_rmatvec(sens)[source]

Adjoint filter for sensitivity backpropagation

Notes

  • Uses cone filter: w(d) = (r_min - d) / r_min

  • Handles non-square elements (dx != dy) automatically

  • Typical r_min values: 1.5-3.0 element units

  • Larger r_min produces smoother, coarser designs

Examples

>>> from pyFANTOM.CPU import StructuredMesh2D, StructuredFilter2D
>>> mesh = StructuredMesh2D(nx=128, ny=64, lx=2.0, ly=1.0)
>>> filter = StructuredFilter2D(mesh=mesh, r_min=2.0)
>>> rho_smooth = filter.dot(rho)
__init__(mesh: StructuredMesh2D, r_min)[source]
class pyFANTOM.geom.CPU._filters.StructuredFilter3D(mesh: StructuredMesh3D, r_min)[source]

Bases: FilterKernel

Convolution-based density filter for 3D structured meshes.

Implements smoothing filter for topology optimization to ensure manufacturability and mesh-independent solutions. Uses efficient convolution with precomputed weights and offsets for structured grids.

Parameters:

  • mesh (StructuredMesh3D) – 3D structured mesh

  • r_min (float) – Filter radius in element units (e.g., r_min=1.5 includes elements within 1.5 element widths)

shape

Filter matrix dimensions (n_elements, n_elements)

Type:

tuple

weights

Precomputed filter weights for neighbor elements

Type:

ndarray

offsets

Element index offsets for neighbors within filter radius

Type:

ndarray

normalizer

Per-element normalization factors for adjoint operation

Type:

ndarray

dot(rho)

Apply forward filter: filtered_rho = H @ rho

_rmatvec(sens)[source]

Apply adjoint filter for sensitivity backpropagation: H^T @ sens

Notes

  • Filter uses cone-shaped kernel: w(d) = (r_min - d) / r_min for d < r_min

  • Automatically handles anisotropic elements (dx != dy != dz)

  • Adjoint includes proper normalization for optimization gradients

  • Memory-efficient: stores only unique weights and offsets, not full matrix

Examples

>>> from pyFANTOM.CPU import StructuredMesh3D, StructuredFilter3D
>>> mesh = StructuredMesh3D(nx=32, ny=32, nz=16, lx=1.0, ly=1.0, lz=0.5)
>>> filter = StructuredFilter3D(mesh=mesh, r_min=1.5)
>>> rho_filtered = filter.dot(rho_raw)
__init__(mesh: StructuredMesh3D, r_min)[source]
class pyFANTOM.geom.CPU._filters.GeneralFilter(mesh: GeneralMesh, r_min)[source]

Bases: FilterKernel

KDTree-based density filter for unstructured meshes.

Implements density filtering for general unstructured meshes using spatial search to identify neighbors. Stores explicit filter matrix (sparse) unlike structured filters.

Parameters:

  • mesh (GeneralMesh) – General unstructured mesh

  • r_min (float) – Filter radius in physical units (not element units)

kernel

Explicit filter matrix, shape (n_elements, n_elements)

Type:

ndarray

shape

Filter matrix dimensions

Type:

tuple

dot(rho)

Forward filter: H @ rho

_rmatvec(sens)[source]

Adjoint filter: H^T @ sens

Notes

  • Uses spatial search (KDTree internally) to build filter

  • Stores explicit sparse matrix (memory-intensive for large meshes)

  • r_min is in physical units, not element-relative

  • Automatically detects 2D vs 3D from mesh

  • Slower initialization than StructuredFilter but handles arbitrary topologies

Examples

>>> from pyFANTOM.CPU import GeneralMesh, GeneralFilter
>>> mesh = GeneralMesh(nodes, elements)  # Triangular mesh
>>> filter = GeneralFilter(mesh=mesh, r_min=0.05)  # Physical radius
>>> rho_filtered = filter.dot(rho)
__init__(mesh: GeneralMesh, r_min)[source]
class pyFANTOM.FiniteElement.CPU.FiniteElement.FiniteElement(mesh: StructuredMesh2D | StructuredMesh3D | GeneralMesh, kernel: StructuredStiffnessKernel | GeneralStiffnessKernel | UniformStiffnessKernel, solver: CHOLMOD | CG | BiCGSTAB | GMRES | SPLU | SPSOLVE | MultiGrid)[source]

Bases: FiniteElement

Finite element analysis engine managing boundary conditions, forces, and solution.

Central class coordinating mesh, stiffness assembly kernel, and linear solver for FEA. Provides convenient methods for applying boundary conditions, forces, solving, and visualization.

Parameters:
mesh

Associated mesh

Type:

Mesh

kernel

Stiffness assembly kernel

Type:

StiffnessKernel

solver

Linear solver

Type:

Solver

rhs

Right-hand side force vector, shape (n_nodes * dof,)

Type:

ndarray

dof

Degrees of freedom per node (2 for 2D, 3 for 3D)

Type:

int

is_3D

True for 3D problems

Type:

bool

add_dirichlet_boundary_condition(node_ids=None, positions=None, dofs=None, rhs=0)[source]

Apply fixed displacement boundary conditions

add_point_forces(forces, node_ids=None, positions=None)[source]

Apply point loads

reset_forces()[source]

Clear all forces

reset_dirichlet_boundary_conditions()[source]

Remove all boundary conditions

solve(rho=None)[source]

Solve linear system K(rho) @ U = F, returns (U, residual)

visualize_problem(**kwargs)[source]

Plot mesh with BCs and loads

visualize_density(rho, \*\*kwargs)[source]

Plot optimization result

visualize_field(field, \*\*kwargs)[source]

Plot displacement or stress field

Notes

  • Boundary conditions modify kernel.constraints boolean array

  • Can specify nodes by index (node_ids) or coordinates (positions with KDTree search)

  • Dirichlet BC handling: rows/columns zeroed, diagonal set to 1

  • solve() returns residual: ||K@U - F|| / ||F|| for validation

Examples

>>> from pyFANTOM.CPU import *
>>> mesh = StructuredMesh2D(nx=64, ny=32, lx=2.0, ly=1.0)
>>> kernel = StructuredStiffnessKernel(mesh=mesh)
>>> solver = CHOLMOD(kernel=kernel)
>>> FE = FiniteElement(mesh=mesh, kernel=kernel, solver=solver)
>>>
>>> # Fix left edge
>>> left_nodes = np.where(mesh.nodes[:, 0] < 1e-6)[0]
>>> FE.add_dirichlet_boundary_condition(node_ids=left_nodes, dofs=np.array([[1,1]]), rhs=0)
>>>
>>> # Apply downward load on right edge
>>> right_nodes = np.where(np.abs(mesh.nodes[:, 0] - 2.0) < 1e-6)[0]
>>> FE.add_point_forces(node_ids=right_nodes, forces=np.array([[0, -1.0]]))
>>>
>>> # Solve
>>> U, residual = FE.solve(rho=np.ones(len(mesh.elements)) * 0.5)
>>> print(f"Residual: {residual:.2e}")
__init__(mesh: StructuredMesh2D | StructuredMesh3D | GeneralMesh, kernel: StructuredStiffnessKernel | GeneralStiffnessKernel | UniformStiffnessKernel, solver: CHOLMOD | CG | BiCGSTAB | GMRES | SPLU | SPSOLVE | MultiGrid)[source]

Initialize finite-element wrapper state.

Subclasses may extend the initializer to accept mesh, kernel or solver objects.

add_dirichlet_boundary_condition(node_ids: ndarray | None = None, positions: ndarray | None = None, dofs: ndarray | None = None, rhs: float | ndarray = 0)[source]

Apply Dirichlet (fixed displacement) boundary conditions.

Parameters:

  • node_ids (ndarray, optional) – Node indices to constrain, shape (n_nodes,)

  • positions (ndarray, optional) – Physical coordinates to constrain (uses KDTree search), shape (n_nodes, spatial_dim)

  • dofs (ndarray, optional) – DOF mask per node, shape (n_nodes, dof) with 1=constrained, 0=free. If None, all DOFs at specified nodes are constrained

  • rhs (float or ndarray, optional) – Prescribed displacement values. Scalar for zero displacement, or array shape (n_nodes, dof)

Notes

  • Provide either node_ids OR positions, not both

  • dofs array: [[1,0]] constrains only x-displacement in 2D

  • Multiple calls accumulate constraints

  • Modifies kernel.constraints boolean array

Examples

>>> # Fix all DOFs at left edge
>>> left_nodes = np.where(mesh.nodes[:, 0] < 1e-6)[0]
>>> FE.add_dirichlet_boundary_condition(node_ids=left_nodes, rhs=0)
>>>
>>> # Fix only y-displacement at bottom
>>> bottom_nodes = np.where(mesh.nodes[:, 1] < 1e-6)[0]
>>> FE.add_dirichlet_boundary_condition(node_ids=bottom_nodes, dofs=np.array([[0,1]]), rhs=0)
>>>
>>> # Prescribe non-zero displacement
>>> FE.add_dirichlet_boundary_condition(node_ids=[10], dofs=np.array([[1,0]]), rhs=np.array([[0.1, 0]]))
add_neumann_boundary_condition(**kwargs)[source]

Add Neumann (natural) boundary condition.

Parameters - condition: backend-specific description of traction/flux loads.

add_point_forces(forces: ndarray, node_ids: ndarray | None = None, positions: ndarray | None = None)[source]

Apply point loads to specified nodes.

Parameters:

  • forces (ndarray) – Force vectors, shape (n_forces, dof) where dof is 2 for 2D or 3 for 3D. For 2D: [[fx, fy]], for 3D: [[fx, fy, fz]]

  • node_ids (ndarray, optional) – Node indices for force application, shape (n_forces,)

  • positions (ndarray, optional) – Physical coordinates for force application (uses KDTree search)

Notes

  • Provide either node_ids OR positions

  • Multiple calls accumulate forces

  • Forces are automatically set to prescribed values at Dirichlet nodes

  • Units should match physics model (e.g., Newtons for SI units)

Examples

>>> # Apply downward force at right edge
>>> right_nodes = np.where(np.abs(mesh.nodes[:, 0] - 2.0) < 1e-6)[0]
>>> FE.add_point_forces(node_ids=right_nodes, forces=np.array([[0, -1.0]]))
>>>
>>> # Apply force at specific location
>>> FE.add_point_forces(positions=np.array([[1.0, 0.5]]), forces=np.array([[10.0, 0]]))
reset_forces()[source]

Clear all applied forces.

Sets the right-hand side force vector to zero. Useful when reconfiguring loading conditions or starting a new analysis.

reset_dirichlet_boundary_conditions()[source]

Remove all Dirichlet boundary conditions.

Clears all fixed displacement constraints and resets the kernel’s constraint state. Useful when reconfiguring boundary conditions.

visualize_problem(ax=None, **kwargs)[source]

Visualize mesh with boundary conditions and loads.

Parameters:

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • **kwargs – Additional arguments passed to plot_problem_2D/3D

Returns:

Plot object (2D: matplotlib axes, 3D: k3d plot)

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Shows: - Mesh elements - Fixed nodes (Dirichlet BCs) as markers - Applied forces as arrows/glyphs - For 3D: Interactive k3d plot - For 2D: Static matplotlib plot

visualize_density(rho, ax=None, **kwargs)[source]

Visualize density distribution (optimization result).

Parameters:

  • rho (ndarray) – Element densities, shape (n_elements,) or (n_elements, n_materials)

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • **kwargs – Additional arguments passed to plot_problem_2D/3D

Returns:

Plot object (2D: matplotlib axes, 3D: k3d plot)

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Color-codes elements by density value (0=void, 1=solid). For multi-material problems, pass rho with shape (n_elements, n_materials).

visualize_field(field, ax=None, rho=None, **kwargs)[source]

Visualize scalar or vector field (displacement, stress, strain, etc.).

Parameters:

  • field (ndarray) – Field values to plot. Shape depends on field type: - Scalar: (n_elements,) or (n_nodes,) - Vector: (n_nodes, dof) for displacement-like fields

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • rho (ndarray, optional) – Density mask to hide void regions, shape (n_elements,)

  • **kwargs – Additional arguments passed to plot_field_2D/3D

Returns:

Plot object (2D: matplotlib axes, 3D: k3d plot)

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Common fields: - Displacement: U.reshape(-1, dof) - Von Mises stress: from Problem.FEA() output - Strain energy: from Problem.FEA() output - If rho provided, elements with rho < 0.5 are hidden

solve(rho=None)[source]

Solve finite element system: K(rho) @ U = F.

Parameters:

rho (ndarray, optional) – Element density variables, shape (n_elements,). If None, uses rho=1 (full density)

Returns:

  • U (ndarray) – Displacement vector, shape (n_nodes * dof,). Interleaved DOFs: [ux0, uy0, ux1, uy1, …]

  • residual (float) – Normalized residual ||K@U - F|| / ||F|| for solution validation

Notes

  • residual < 1e-5 indicates accurate solution

  • residual > 1e-2 indicates ill-conditioned system (check Problem.ill_conditioned())

  • Solver reuses factorization if available (e.g., CHOLMOD)

  • For topology optimization, rho represents material distribution

Examples

>>> rho = np.ones(len(mesh.elements)) * 0.5  # 50% density everywhere
>>> U, residual = FE.solve(rho=rho)
>>> if residual > 1e-3:
>>>     print("Warning: Poor convergence!")
>>>
>>> # Extract displacements
>>> U_xy = U.reshape(-1, 2)  # Shape: (n_nodes, 2) for 2D
>>> u_x = U_xy[:, 0]
>>> u_y = U_xy[:, 1]
class pyFANTOM.FiniteElement.CPU.NLFiniteElement.NLFiniteElement(mesh: StructuredMesh2D | StructuredMesh3D | GeneralMesh, kernel: StructuredStiffnessKernel | GeneralStiffnessKernel | UniformStiffnessKernel, solver: CHOLMOD | CG | BiCGSTAB | GMRES | SPLU | SPSOLVE | MultiGrid, physics: NLElasticity)[source]

Bases: FiniteElement

Nonlinear finite element analysis engine for geometrically nonlinear problems.

Extends FiniteElement with Newton-Raphson solver for large deformation analysis. Handles geometrically nonlinear elasticity where element stiffness matrices depend on current deformation state. Requires NLElasticity physics and NLUniformStiffnessKernel.

Parameters:
mesh

Associated mesh

Type:

Mesh

kernel

Nonlinear stiffness assembly kernel

Type:

NLUniformStiffnessKernel

solver

Linear solver

Type:

Solver

rhs

Right-hand side force vector, shape (n_nodes * dof,)

Type:

ndarray

dof

Degrees of freedom per node (2 for 2D, 3 for 3D)

Type:

int

is_3D

True for 3D problems

Type:

bool

state

Nonlinear state manager tracking deformation and internal forces

Type:

NLState

lagrangeMult

Lagrange multipliers from adjoint solve (for sensitivity analysis)

Type:

ndarray

last_dR_Drho

Derivative of residual w.r.t. density (for sensitivity analysis)

Type:

ndarray

add_dirichlet_boundary_condition(node_ids=None, positions=None, dofs=None, rhs=0)[source]

Apply fixed displacement boundary conditions

add_point_forces(forces, node_ids=None, positions=None)[source]

Apply point loads

reset_forces()[source]

Clear all forces

reset_dirichlet_boundary_conditions()[source]

Remove all boundary conditions

solveNL(rho=None)[source]

Solve nonlinear system using Newton-Raphson, returns (U, residual)

visualize_problem(**kwargs)[source]

Plot mesh with BCs and loads

visualize_density(rho, \*\*kwargs)[source]

Plot optimization result

visualize_deformed_mesh(rho=None, \*\*kwargs)[source]

Plot deformed mesh configuration

visualize_field(field, \*\*kwargs)[source]

Plot displacement or stress field

Notes

Nonlinear Analysis: - Uses Newton-Raphson method with load stepping - Tangent stiffness matrix updated each iteration - Load stepping: 4 load steps by default for convergence - Convergence based on out-of-balance work

Sensitivity Analysis: - Solves adjoint system after equilibrium - Computes dR/drho for gradient computation - More expensive than linear analysis

Use Cases: - Large deformation problems - Buckling-sensitive structures - Compliant mechanisms - Soft robotics

Examples

>>> from pyFANTOM.CPU import *
>>> from pyFANTOM import NLElasticity
>>> # Setup nonlinear FEA
>>> physics = NLElasticity(E=1.0, nu=0.3)
>>> mesh = StructuredMesh2D(nx=64, ny=32, lx=2.0, ly=1.0, physics=physics)
>>> kernel = NLUniformStiffnessKernel(mesh=mesh)
>>> solver = CG(kernel=kernel)
>>> FE = NLFiniteElement(mesh=mesh, kernel=kernel, solver=solver, physics=physics)
>>>
>>> # Apply BCs and loads
>>> FE.add_dirichlet_boundary_condition(node_ids=fixed_nodes, rhs=0)
>>> FE.add_point_forces(node_ids=load_nodes, forces=np.array([[0, -1.0]]))
>>>
>>> # Solve nonlinear system
>>> U, residual = FE.solveNL(rho=np.ones(len(mesh.elements)) * 0.5)
>>> print(f"Residual: {residual:.2e}")
__init__(mesh: StructuredMesh2D | StructuredMesh3D | GeneralMesh, kernel: StructuredStiffnessKernel | GeneralStiffnessKernel | UniformStiffnessKernel, solver: CHOLMOD | CG | BiCGSTAB | GMRES | SPLU | SPSOLVE | MultiGrid, physics: NLElasticity)[source]

Initialize finite-element wrapper state.

Subclasses may extend the initializer to accept mesh, kernel or solver objects.

add_dirichlet_boundary_condition(node_ids: ndarray | None = None, positions: ndarray | None = None, dofs: ndarray | None = None, rhs: float | ndarray = 0)[source]

Apply Dirichlet (fixed displacement) boundary conditions.

Parameters:

  • node_ids (ndarray, optional) – Node indices to constrain, shape (n_nodes,)

  • positions (ndarray, optional) – Physical coordinates to constrain (uses KDTree search), shape (n_nodes, spatial_dim)

  • dofs (ndarray, optional) – DOF mask per node, shape (n_nodes, dof) with 1=constrained, 0=free. If None, all DOFs at specified nodes are constrained

  • rhs (float or ndarray, optional) – Prescribed displacement values. Scalar for zero displacement, or array shape (n_nodes, dof)

Notes

  • Provide either node_ids OR positions, not both

  • dofs array: [[1,0]] constrains only x-displacement in 2D

  • Multiple calls accumulate constraints

  • Modifies kernel.constraints boolean array

  • Same interface as linear FiniteElement

add_neumann_boundary_condition(**kwargs)[source]

Apply Neumann boundary conditions (not implemented).

Raises:

NotImplementedError – Always raised. Neumann boundary conditions not yet implemented for nonlinear FEA.

add_point_forces(forces: ndarray, node_ids: ndarray | None = None, positions: ndarray | None = None)[source]

Apply point loads to specified nodes.

Parameters:

  • forces (ndarray) – Force vectors, shape (n_forces, dof) where dof is 2 for 2D or 3 for 3D. For 2D: [[fx, fy]], for 3D: [[fx, fy, fz]]

  • node_ids (ndarray, optional) – Node indices for force application, shape (n_forces,)

  • positions (ndarray, optional) – Physical coordinates for force application (uses KDTree search)

Notes

  • Provide either node_ids OR positions

  • Multiple calls accumulate forces

  • Forces are automatically set to prescribed values at Dirichlet nodes

  • Units should match physics model (e.g., Newtons for SI units)

  • Same interface as linear FiniteElement

reset_forces()[source]

Clear all applied forces.

Sets the right-hand side force vector to zero. Useful when reconfiguring loading conditions or starting a new analysis.

reset_dirichlet_boundary_conditions()[source]

Remove all Dirichlet boundary conditions.

Clears all fixed displacement constraints and resets the kernel’s constraint state. Useful when reconfiguring boundary conditions.

visualize_problem(ax=None, **kwargs)[source]

Visualize mesh with boundary conditions and loads.

Parameters:

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • **kwargs – Additional arguments passed to plot_problem_2D/3D

Returns:

Plot object (2D: matplotlib axes, 3D: k3d plot)

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Shows: - Mesh elements - Fixed nodes (Dirichlet BCs) as markers - Applied forces as arrows/glyphs - For 3D: Interactive k3d plot - For 2D: Static matplotlib plot

visualize_solution(ax=None, **kwargs)[source]

Visualize solution (mesh with BCs and loads).

Parameters:

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • **kwargs – Additional arguments passed to plot_problem_2D

Returns:

Plot object (2D only, 3D raises NotImplementedError)

Return type:

matplotlib.axes.Axes

Notes

Currently only supports 2D visualization. For 3D, use visualize_problem().

visualize_density(rho, ax=None, **kwargs)[source]

Visualize density distribution (optimization result).

Parameters:

  • rho (ndarray) – Element densities, shape (n_elements,) or (n_elements, n_materials)

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • **kwargs – Additional arguments passed to plot_problem_2D/3D

Returns:

Plot object (2D: matplotlib axes, 3D: k3d plot)

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Color-codes elements by density value (0=void, 1=solid). For multi-material problems, pass rho with shape (n_elements, n_materials).

visualize_deformed_mesh(rho=None, ax=None, **kwargs)[source]

Visualize deformed mesh configuration.

Parameters:

  • rho (ndarray, optional) – Element densities for masking void regions, shape (n_elements,). If None, uses full density (rho=1.0)

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • **kwargs – Additional arguments passed to plot_problem_2D_wVoid

Returns:

Plot object (2D only, 3D raises NotImplementedError)

Return type:

matplotlib.axes.Axes

Notes

Shows mesh in deformed configuration (nodes + displacements from solveNL()). Useful for visualizing large deformations in nonlinear analysis. Currently only supports 2D visualization.

visualize_field(field, ax=None, rho=None, **kwargs)[source]

Visualize scalar or vector field (displacement, stress, strain, etc.).

Parameters:

  • field (ndarray) – Field values to plot. Shape depends on field type: - Scalar: (n_elements,) or (n_nodes,) - Vector: (n_nodes, dof) for displacement-like fields

  • ax (matplotlib.axes.Axes, optional) – Existing axes to plot on (2D only). If None, creates new figure

  • rho (ndarray, optional) – Density mask to hide void regions, shape (n_elements,)

  • **kwargs – Additional arguments passed to plot_field_2D/3D

Returns:

Plot object (2D: matplotlib axes, 3D: k3d plot)

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Common fields: - Displacement: U.reshape(-1, dof) - Von Mises stress: from Problem.FEA() output - Strain energy: from Problem.FEA() output - If rho provided, elements with rho < 0.5 are hidden

solve(rho=None)[source]

Solve linear system (for compatibility, calls solveNL).

Parameters:

rho (ndarray, optional) – Element density variables, shape (n_elements,). If None, uses rho=1 (full density)

Returns:

  • U (ndarray) – Displacement vector, shape (n_nodes * dof,)

  • residual (float) – Normalized residual ||K@U - F|| / ||F||

Notes

For compatibility with linear FiniteElement interface. In practice, use solveNL() for nonlinear analysis. This method may not work correctly for large deformations.

solveNL(rho=None)[source]

Solve nonlinear system using Newton-Raphson method with load stepping.

Parameters:

rho (ndarray, optional) – Element density variables, shape (n_elements,). If None, uses rho=1 (full density)

Returns:

  • U (ndarray) – Total displacement vector after convergence, shape (n_nodes * dof,)

  • residual (float) – Final normalized residual ||K@U - F|| / ||F||

Notes

Newton-Raphson Method: - Uses load stepping: 4 load steps by default - Each load step applies incremental force: F_step = (step+1)/nLoadSteps * F_total - Up to 50 Newton iterations per load step - Convergence: |OoBWork_current / OoBWork_initial| < 5e-5

Out-of-Balance Work: - OoBWork = K@U - F (residual force vector) - Measures equilibrium error - Constrained DOFs excluded from convergence check

Tangent Stiffness: - Updated each Newton iteration via state._KTan() - Computed from current deformation state - Passed to kernel via set_Ktan()

Adjoint Preparation: - After equilibrium, solves for Lagrange multipliers - Stores dR/drho for sensitivity analysis - Required for gradient computation in optimization

State Management: - Updates state.tDeltatUk (incremental displacement) - Updates state.tU_global (total displacement) - Tracks element stress and deformation

Raises:

ValueError – If rho has wrong shape or solver doesn’t converge within max iterations

Examples

>>> U, residual = FE.solveNL(rho=np.ones(nel) * 0.5)
>>> print(f"Converged with residual: {residual:.2e}")
>>> # Access final displacement
>>> U_reshaped = U.reshape(-1, 2)  # For 2D: (n_nodes, 2)
class pyFANTOM.stiffness.CPU._FEA.StructuredStiffnessKernel(mesh: StructuredMesh)[source]

Bases: StiffnessKernel

Optimized stiffness matrix assembly for structured meshes.

Uses node-basis assembly with a single precomputed element stiffness matrix for maximum efficiency on uniform grid meshes. Supports matrix-free operations (matvec) and explicit CSR construction.

Parameters:

mesh (StructuredMesh) – StructuredMesh2D or StructuredMesh3D with uniform elements

K_single

Single element stiffness matrix used for all elements

Type:

ndarray

shape

Global stiffness matrix dimensions (n_dof, n_dof)

Type:

tuple

constraints

Boolean array marking constrained DOFs

Type:

ndarray (bool)

has_cons

True if boundary conditions have been applied

Type:

bool

non_con_map

Index map to non-constrained DOFs for reduced system

Type:

ndarray

set_rho(rho)

Set design variables for subsequent matrix-free operations

dot(rhs)[source]

Matrix-vector product K @ rhs (requires rho to be set)

construct(rho)[source]

Build explicit CSR matrix representation

diagonal(rho=None)[source]

Extract diagonal of stiffness matrix

add_constraints(dof_indices)[source]

Apply Dirichlet boundary conditions

process_grad(U)[source]

Compute element-wise sensitivities dK/drho : U for adjoint method

Notes

  • Fastest kernel for structured meshes (2-3x faster than UniformStiffnessKernel)

  • Automatically handles constrained DOFs by zeroing rows/columns

  • Matrix-free dot() operation uses optimized numba kernels

  • Use construct() only when explicit matrix is needed (e.g., for direct solvers)

Examples

>>> from pyFANTOM.CPU import StructuredMesh2D, StructuredStiffnessKernel
>>> mesh = StructuredMesh2D(nx=64, ny=32, lx=2.0, ly=1.0)
>>> kernel = StructuredStiffnessKernel(mesh=mesh)
>>> kernel.add_constraints(np.array([0, 1, 2]))  # Fix first node
>>> rho = np.ones(len(mesh.elements)) * 0.5
>>> kernel.set_rho(rho)
>>> u = kernel.dot(f)  # Matrix-free K @ f
__init__(mesh: StructuredMesh)[source]
diagonal(rho=None)[source]

Extract diagonal of stiffness matrix.

Parameters:

rho (ndarray, optional) – Design variables. If None, uses self.rho (must be set via set_rho())

Returns:

Diagonal entries, shape (n_dof,)

Return type:

ndarray

Raises:

ValueError – If rho is None and has_rho is False

set_constraints(constraints)[source]

Set Dirichlet boundary conditions (replaces existing constraints).

Parameters:

constraints (ndarray) – DOF indices to constrain, shape (n_constrained_dofs,)

Notes

Clears all existing constraints and sets only the specified DOFs. Use add_constraints() to add constraints incrementally.

add_constraints(constraints)[source]

Add Dirichlet boundary conditions (accumulates with existing constraints).

Parameters:

constraints (ndarray) – DOF indices to constrain, shape (n_constrained_dofs,)

Notes

Adds to existing constraints. Use set_constraints() to replace all constraints.

process_grad(U)[source]

Compute element-wise sensitivities for adjoint method.

Parameters:

U (ndarray) – Displacement solution vector, shape (n_dof,)

Returns:

Element-wise sensitivities dK/drho : U, shape (n_elements,)

Return type:

ndarray

Notes

Computes -U^T @ dK/drho @ U for each element, used in compliance sensitivity analysis. Result is per-element contribution to objective gradient.

construct(rho)[source]

Build explicit CSR sparse matrix representation.

Parameters:

rho (ndarray) – Design variables, shape (n_elements,)

Returns:

Sparse stiffness matrix in CSR format

Return type:

csr_matrix

Notes

  • Memory-intensive: O(nnz) storage

  • Required for direct solvers (CHOLMOD, SPLU)

  • Reuses sparsity pattern on subsequent calls (faster)

dot(rhs)[source]

Matrix-vector or matrix-matrix product K(rho) @ rhs.

Parameters:

rhs (ndarray or csr_matrix) – Right-hand side vector or matrix - Vector: shape (n_dof,) - Matrix: shape (n_dof, n_cols)

Returns:

Result of K @ rhs, same shape as rhs

Return type:

ndarray or csr_matrix

Raises:

  • ValueError – If rho has not been set (call set_rho() first)

  • ValueError – If rhs shape doesn’t match matrix dimensions

Notes

Requires rho to be set via set_rho() or passed to construct().

class pyFANTOM.stiffness.CPU._FEA.GeneralStiffnessKernel(mesh: GeneralMesh)[source]

Bases: StiffnessKernel

Stiffness matrix assembly for fully heterogeneous meshes.

Most general kernel supporting meshes with elements of varying sizes and topologies. Uses flat storage with pointer arrays for memory-efficient representation.

Parameters:

mesh (GeneralMesh) – GeneralMesh with potentially heterogeneous elements (mesh.is_uniform can be False)

K_flat

Flattened element stiffness matrices

Type:

ndarray

K_ptr

Pointer array for accessing K_flat: K_e = K_flat[K_ptr[e]:K_ptr[e+1]]

Type:

ndarray

elements_flat

Flattened element connectivity

Type:

ndarray

elements_ptr

Pointer array for accessing elements

Type:

ndarray

element_sizes

Number of nodes per element, shape (n_elements,)

Type:

ndarray

Same interface as StructuredStiffnessKernel

Notes

  • Supports mixed element types (triangles + quads, tets + hexes, etc.)

  • Flat storage with pointers enables variable element sizes

  • Slower than Uniform/Structured kernels due to indirection

  • Automatically used if mesh.is_uniform == False

Examples

>>> from pyFANTOM.CPU import GeneralMesh, GeneralStiffnessKernel
>>> # Mixed triangular and quad mesh
>>> nodes = np.array([[0,0], [1,0], [0,1], [1,1], [2,0]])
>>> elements = [np.array([0,1,2]), np.array([1,4,3,2])]  # Triangle + Quad
>>> mesh = GeneralMesh(nodes, elements)
>>> kernel = GeneralStiffnessKernel(mesh=mesh)
__init__(mesh: GeneralMesh)[source]
diagonal(rho=None)[source]

Extract diagonal of stiffness matrix.

Parameters:

rho (ndarray, optional) – Design variables. If None, uses self.rho (must be set via set_rho())

Returns:

Diagonal entries, shape (n_dof,)

Return type:

ndarray

Raises:

ValueError – If rho is None and has_rho is False

set_constraints(constraints)[source]

Set Dirichlet boundary conditions (replaces existing constraints).

Parameters:

constraints (ndarray) – DOF indices to constrain, shape (n_constrained_dofs,)

Notes

Clears all existing constraints and sets only the specified DOFs. Use add_constraints() to add constraints incrementally.

add_constraints(constraints)[source]

Add Dirichlet boundary conditions (accumulates with existing constraints).

Parameters:

constraints (ndarray) – DOF indices to constrain, shape (n_constrained_dofs,)

Notes

Adds to existing constraints. Use set_constraints() to replace all constraints.

process_grad(U)[source]

Compute element-wise sensitivities for adjoint method.

Parameters:

U (ndarray) – Displacement solution vector, shape (n_dof,)

Returns:

Element-wise sensitivities dK/drho : U, shape (n_elements,)

Return type:

ndarray

Notes

Computes -U^T @ dK/drho @ U for each element, used in compliance sensitivity analysis. Result is per-element contribution to objective gradient.

construct(rho)[source]

Build explicit CSR sparse matrix representation.

Parameters:

rho (ndarray) – Design variables, shape (n_elements,)

Returns:

Sparse stiffness matrix in CSR format

Return type:

csr_matrix

Notes

  • Memory-intensive: O(nnz) storage

  • Required for direct solvers (CHOLMOD, SPLU)

  • Reuses sparsity pattern on subsequent calls (faster)

dot(rhs)[source]

Matrix-vector or matrix-matrix product K(rho) @ rhs.

Parameters:

rhs (ndarray or csr_matrix) – Right-hand side vector or matrix - Vector: shape (n_dof,) - Matrix: shape (n_dof, n_cols)

Returns:

Result of K @ rhs, same shape as rhs

Return type:

ndarray or csr_matrix

Raises:

  • ValueError – If rho has not been set (call set_rho() first)

  • ValueError – If rhs shape doesn’t match matrix dimensions

Notes

Requires rho to be set via set_rho() or passed to construct().

class pyFANTOM.stiffness.CPU._FEA.UniformStiffnessKernel(mesh: GeneralMesh)[source]

Bases: StiffnessKernel

Stiffness matrix assembly for general uniform meshes.

Handles unstructured meshes where all elements have the same number of nodes and identical geometry, using element-specific stiffness matrices. More general than StructuredStiffnessKernel but less efficient.

Parameters:

mesh (GeneralMesh) – GeneralMesh with uniform element topology (mesh.is_uniform must be True)

Ks

Element stiffness matrices, shape (n_elements, dof_per_element, dof_per_element)

Type:

ndarray

shape

Global stiffness matrix dimensions

Type:

tuple

constraints

Boolean array for constrained DOFs

Type:

ndarray (bool)

Same interface as StructuredStiffnessKernel

Notes

  • Requires mesh.is_uniform == True (all elements same size)

  • Stores per-element stiffness matrices (more memory than StructuredStiffnessKernel)

  • Use for unstructured uniform meshes (e.g., triangular or tetrahedral)

  • For structured grids, prefer StructuredStiffnessKernel for better performance

Examples

>>> from pyFANTOM.CPU import GeneralMesh, UniformStiffnessKernel
>>> # Triangular mesh
>>> nodes = np.array([[0,0], [1,0], [0,1], [1,1]])
>>> elements = np.array([[0,1,2], [1,3,2]])
>>> mesh = GeneralMesh(nodes, elements)
>>> kernel = UniformStiffnessKernel(mesh=mesh)
__init__(mesh: GeneralMesh)[source]
diagonal(rho=None)[source]

Extract diagonal of stiffness matrix.

Parameters:

rho (ndarray, optional) – Design variables. If None, uses self.rho (must be set via set_rho())

Returns:

Diagonal entries, shape (n_dof,)

Return type:

ndarray

Raises:

ValueError – If rho is None and has_rho is False

set_constraints(constraints)[source]

Set Dirichlet boundary conditions (replaces existing constraints).

Parameters:

constraints (ndarray) – DOF indices to constrain, shape (n_constrained_dofs,)

Notes

Clears all existing constraints and sets only the specified DOFs. Use add_constraints() to add constraints incrementally.

add_constraints(constraints)[source]

Add Dirichlet boundary conditions (accumulates with existing constraints).

Parameters:

constraints (ndarray) – DOF indices to constrain, shape (n_constrained_dofs,)

Notes

Adds to existing constraints. Use set_constraints() to replace all constraints.

process_grad(U)[source]

Compute element-wise sensitivities for adjoint method.

Parameters:

U (ndarray) – Displacement solution vector, shape (n_dof,)

Returns:

Element-wise sensitivities dK/drho : U, shape (n_elements,)

Return type:

ndarray

Notes

Computes -U^T @ dK/drho @ U for each element, used in compliance sensitivity analysis. Result is per-element contribution to objective gradient.

construct(rho)[source]

Build explicit CSR sparse matrix representation.

Parameters:

rho (ndarray) – Design variables, shape (n_elements,)

Returns:

Sparse stiffness matrix in CSR format

Return type:

csr_matrix

Notes

  • Memory-intensive: O(nnz) storage

  • Required for direct solvers (CHOLMOD, SPLU)

  • Reuses sparsity pattern on subsequent calls (faster)

dot(rhs)[source]

Matrix-vector or matrix-matrix product K(rho) @ rhs.

Parameters:

rhs (ndarray or csr_matrix) – Right-hand side vector or matrix - Vector: shape (n_dof,) - Matrix: shape (n_dof, n_cols)

Returns:

Result of K @ rhs, same shape as rhs

Return type:

ndarray or csr_matrix

Raises:

  • ValueError – If rho has not been set (call set_rho() first)

  • ValueError – If rhs shape doesn’t match matrix dimensions

Notes

Requires rho to be set via set_rho() or passed to construct().

class pyFANTOM.solvers.CPU._solvers.CHOLMOD(kernel: StiffnessKernel)[source]

Bases: Solver

Direct sparse Cholesky solver using CHOLMOD (scikit-sparse).

Most efficient direct solver for symmetric positive definite systems. Uses supernodal Cholesky factorization with reusable symbolic factorization. Recommended for problems with <3M DOF on CPU.

Parameters:

kernel (StiffnessKernel) – Stiffness assembly kernel (must have shape[0] <= 3M DOF)

factor

CHOLMOD factorization object (reused across solves)

Type:

cholmod.Factor

factorized

True after first factorization

Type:

bool

solve(rhs, rho=None, \*\*kwargs)[source]

Solve K(rho) @ U = rhs using Cholesky factorization

reset()[source]

Clear factorization to force reinitialization

initialize()[source]

Compute symbolic factorization

Notes

  • Fastest solver for small-medium problems (<500k DOF)

  • Requires MKL-optimized scikit-sparse for best performance

  • Reuses symbolic factorization across iterations (only numerical refactorization)

  • Memory intensive: ~O(n^1.5) for 2D, ~O(n^2) for 3D

  • For >3M DOF, use iterative solvers (CG, MultiGrid)

Examples

>>> from pyFANTOM.CPU import StructuredStiffnessKernel, CHOLMOD
>>> solver = CHOLMOD(kernel=kernel)
>>> U, residual = solver.solve(rhs=F, rho=rho)
>>> print(f"Residual: {residual:.2e}")
__init__(kernel: StiffnessKernel)[source]
reset()[source]

Reset solver state.

Clears internal state such as factorizations, preconditioners, or iteration history. Useful when mesh or boundary conditions change significantly.

initialize()[source]
solve(rhs, rho=None, **kwargs)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.solvers.CPU._solvers.CG(kernel: StiffnessKernel, maxiter=1000, tol=1e-05, matrix_free=None)[source]

Bases: Solver

Conjugate Gradient iterative solver.

Memory-efficient Krylov subspace solver for symmetric positive definite systems. Supports both matrix-free and explicit matrix modes. Best for well-conditioned problems.

Parameters:

  • kernel (StiffnessKernel) – Stiffness assembly kernel

  • maxiter (int, optional) – Maximum iterations (default: 1000)

  • tol (float, optional) – Relative convergence tolerance (default: 1e-5)

  • matrix_free (bool, optional) – Use matrix-free operations (True) or explicit matrix (False). Auto-detected: True for StructuredStiffnessKernel, False for General/Uniform

last_x0

Last solution (used as initial guess for warm-starting)

Type:

ndarray

solve(rhs, rho=None, use_last=True)[source]

Solve K(rho) @ U = rhs iteratively

Notes

  • Matrix-free mode: O(n) memory, slower per iteration

  • Explicit mode: O(nnz) memory, faster per iteration for unstructured meshes

  • Warm-starting (use_last=True) accelerates successive solves

  • Convergence depends on condition number: κ(K)

  • For ill-conditioned problems, consider preconditioners or MultiGrid

Examples

>>> solver = CG(kernel=kernel, maxiter=500, tol=1e-6)
>>> U, residual = solver.solve(rhs=F, rho=rho, use_last=True)
__init__(kernel: StiffnessKernel, maxiter=1000, tol=1e-05, matrix_free=None)[source]
solve(rhs, rho=None, use_last=True)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.solvers.CPU._solvers.BiCGSTAB(kernel: StiffnessKernel, maxiter=1000, tol=1e-05, matrix_free=None)[source]

Bases: Solver

Biconjugate Gradient Stabilized iterative solver.

More robust than CG for non-symmetric or poorly conditioned systems. Rarely needed for standard elasticity problems but useful for coupled physics or unusual BCs.

Parameters:

  • kernel (StiffnessKernel) – Stiffness assembly kernel

  • maxiter (int, optional) – Maximum iterations (default: 1000)

  • tol (float, optional) – Relative convergence tolerance (default: 1e-5)

  • matrix_free (bool, optional) – Use matrix-free operations or explicit matrix

Notes

  • More expensive per iteration than CG (~2x work)

  • Better convergence for indefinite systems

  • Supports warm-starting like CG

  • For elasticity, CG is usually sufficient and faster

Examples

>>> solver = BiCGSTAB(kernel=kernel, tol=1e-5)
>>> U, residual = solver.solve(rhs=F, rho=rho)
__init__(kernel: StiffnessKernel, maxiter=1000, tol=1e-05, matrix_free=None)[source]
solve(rhs, rho=None, use_last=True)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.solvers.CPU._solvers.GMRES(kernel: StiffnessKernel, maxiter=1000, tol=1e-05, matrix_free=None)[source]

Bases: Solver

Generalized Minimal Residual iterative solver.

Most robust Krylov solver for general linear systems. Higher memory usage due to Arnoldi process. Use when CG/BiCGSTAB fail to converge.

Parameters:

  • kernel (StiffnessKernel) – Stiffness assembly kernel

  • maxiter (int, optional) – Maximum iterations (default: 1000)

  • tol (float, optional) – Relative convergence tolerance (default: 1e-5)

  • matrix_free (bool, optional) – Use matrix-free operations or explicit matrix

Notes

  • Memory: O(n * restart), where restart is typically 20

  • More robust than CG for indefinite or non-symmetric systems

  • Slower per iteration than CG

  • For topology optimization, CG is usually preferred

Examples

>>> solver = GMRES(kernel=kernel, maxiter=500)
>>> U, residual = solver.solve(rhs=F, rho=rho)
__init__(kernel: StiffnessKernel, maxiter=1000, tol=1e-05, matrix_free=None)[source]
solve(rhs, rho=None, use_last=True)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.solvers.CPU._solvers.SPLU(kernel: StiffnessKernel)[source]

Bases: Solver

Direct sparse LU solver using SuperLU.

General-purpose direct solver using LU factorization. Works for non-symmetric systems but less efficient than CHOLMOD for SPD matrices. Limited to <3M DOF.

Parameters:

kernel (StiffnessKernel) – Stiffness assembly kernel (must have shape[0] <= 3M DOF)

Notes

  • More general than CHOLMOD (handles non-symmetric matrices)

  • Less efficient than CHOLMOD for elasticity (SPD systems)

  • No factorization reuse (refactorizes every solve)

  • Memory intensive

  • For elasticity, prefer CHOLMOD over SPLU

Examples

>>> solver = SPLU(kernel=kernel)
>>> U, residual = solver.solve(rhs=F, rho=rho)
__init__(kernel: StiffnessKernel)[source]
solve(rhs, rho=None, **kwargs)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.solvers.CPU._solvers.SPSOLVE(kernel: StiffnessKernel)[source]

Bases: Solver

Direct sparse solver using scipy.sparse.linalg.spsolve.

Convenience wrapper around SciPy’s spsolve (auto-selects LU or Cholesky). Simple to use but no factorization reuse. Limited to <3M DOF.

Parameters:

kernel (StiffnessKernel) – Stiffness assembly kernel (must have shape[0] <= 3M DOF)

Notes

  • Auto-detects matrix symmetry and chooses solver

  • No factorization reuse (slow for repeated solves)

  • Simple API but inefficient for optimization

  • For production use, prefer CHOLMOD (better performance)

Examples

>>> solver = SPSOLVE(kernel=kernel)
>>> U, residual = solver.solve(rhs=F, rho=rho)
__init__(kernel: StiffnessKernel)[source]
solve(rhs, rho=None, **kwargs)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.solvers.CPU._solvers.MultiGrid(mesh: StructuredMesh2D | StructuredMesh3D, kernel: StiffnessKernel, maxiter=1000, tol=1e-05, n_smooth=3, omega=0.5, n_level=3, cycle='W', w_level=1, coarse_solver='cholmod', matrix_free=False, low_level_tol=1e-08, low_level_maxiter=5000, min_omega=0.4, omega_boost=1.06)[source]

Bases: Solver

Geometric multigrid solver for structured meshes.

Memory-efficient iterative solver using hierarchical grid refinement. Best choice for large 3D problems (>1M DOF) where direct solvers are impractical. Significantly faster than CG for well-conditioned elasticity problems.

Parameters:

  • mesh (StructuredMesh2D or StructuredMesh3D) – Structured mesh (required for geometric coarsening)

  • kernel (StiffnessKernel) – Stiffness assembly kernel

  • maxiter (int, optional) – Maximum V/W-cycles (default: 1000)

  • tol (float, optional) – Relative convergence tolerance (default: 1e-5)

  • n_smooth (int, optional) – Smoothing iterations per level (default: 3)

  • omega (float, optional) – Damping parameter for Jacobi smoother (default: 0.5)

  • n_level (int, optional) – Number of multigrid levels (default: 3)

  • cycle (str, optional) – Cycle type: ‘V’ or ‘W’ (default: ‘W’)

  • w_level (int or list, optional) – Level(s) for W-cycle recursion (default: 1)

  • coarse_solver (str, optional) – Coarsest level solver: ‘cg’, ‘bicgstab’, ‘gmres’, ‘splu’, ‘spsolve’, ‘cholmod’ (default: ‘splu’)

  • matrix_free (bool, optional) – Use matrix-free operations (default: False, not yet implemented)

  • low_level_tol (float, optional) – Coarse solver tolerance (default: 1e-8)

  • low_level_maxiter (int, optional) – Coarse solver max iterations (default: 5000)

cycle

‘V’ or ‘W’ cycle type

Type:

str

n_level

Number of multigrid levels

Type:

int

PRs

Prolongation/restriction operators

Type:

list

solve(rhs, rho=None, use_last=True)[source]

Solve K(rho) @ U = rhs using multigrid cycles

Notes

  • Best for large 3D problems: Scales O(n) vs O(n^1.5) for CHOLMOD

  • GPU-friendly: Primary solver for CUDA backend

  • W-cycle: More expensive but better convergence than V-cycle

  • Requires structured mesh: Cannot be used with GeneralMesh

  • Typical performance: 5-10x faster than CG for 3D problems >500k DOF

  • Coarse solver: ‘splu’ for small problems, ‘cg’ for very large coarse grids

Examples

>>> from pyFANTOM.CPU import StructuredMesh3D, StructuredStiffnessKernel, MultiGrid
>>> mesh = StructuredMesh3D(nx=64, ny=64, nz=64, lx=1.0, ly=1.0, lz=1.0)
>>> kernel = StructuredStiffnessKernel(mesh=mesh)
>>> solver = MultiGrid(mesh=mesh, kernel=kernel, n_level=4, cycle='W')
>>> U, residual = solver.solve(rhs=F, rho=rho)
>>> print(f"Converged in {solver.maxiter} cycles, residual: {residual:.2e}")
__init__(mesh: StructuredMesh2D | StructuredMesh3D, kernel: StiffnessKernel, maxiter=1000, tol=1e-05, n_smooth=3, omega=0.5, n_level=3, cycle='W', w_level=1, coarse_solver='cholmod', matrix_free=False, low_level_tol=1e-08, low_level_maxiter=5000, min_omega=0.4, omega_boost=1.06)[source]
reset()[source]

Reset solver state.

Clears internal state such as factorizations, preconditioners, or iteration history. Useful when mesh or boundary conditions change significantly.

solve(rhs, rho=None, use_last=True)[source]

Solve linear system K(rho) @ U = rhs.

Parameters:

  • rhs (ndarray) – Right-hand side force vector, shape (n_dof,)

  • rho (ndarray, optional) – Design variables (densities), shape (n_elements,)

  • **kwargs – Additional solver-specific parameters

Returns:

  • U (ndarray) – Solution vector, shape (n_dof,)

  • residual (float) – Relative residual: ||K@U - rhs|| / ||rhs||

Raises:

NotImplementedError – Must be implemented in subclasses

class pyFANTOM.Optimizers.CPU.MMA.MMA(problem: Problem, sub_tol=1e-07, sub_maxiter=100, change_tol=0.0001, fun_tol=1e-06, move=0.5, timer=False)[source]

Bases: Optimizer

Method of Moving Asymptotes (MMA) nonlinear optimizer.

Industry-standard gradient-based optimizer for topology optimization. Handles general nonlinear constraints efficiently via convex approximations. Most robust optimizer in pyFANTOM.

Parameters:

  • problem (Problem) – Optimization problem (e.g., MinimumCompliance)

  • sub_tol (float, optional) – Sub-problem convergence tolerance (default: 1e-7)

  • sub_maxiter (int, optional) – Sub-problem max iterations (default: 100)

  • change_tol (float, optional) – Design variable change convergence tolerance (default: 1e-4)

  • fun_tol (float, optional) – Objective function change tolerance (default: 1e-6)

  • move (float, optional) – Move limit for design variables, range [0,1] (default: 0.5)

  • timer (bool, optional) – Return iteration time if True (default: False)

iteration

Current iteration number

Type:

int

change

L2 norm of design variable change

Type:

float

change_f

Relative objective function change

Type:

float

iter()[source]

Perform one MMA iteration

converged()[source]

Check convergence based on change_tol and fun_tol

logs()[source]

Return dict of optimization metrics

Notes

  • Best for: General constraints, multi-material, multi-constraint problems

  • Convergence: Typically 50-200 iterations for topology optimization

  • Move limit: Smaller = more stable but slower, larger = faster but may oscillate

  • Sub-problem: Convex approximation solved each iteration

  • Dual variables: Automatically updated via Lagrange multipliers

Examples

>>> from pyFANTOM.CPU import MMA, MinimumCompliance
>>> optimizer = MMA(problem=problem, sub_tol=1e-7, fun_tol=1e-6)
>>> for i in range(100):
>>>     optimizer.iter()
>>>     logs = optimizer.logs()
>>>     print(f"Iter {i}: C={logs['objective']:.2e}, Vol={logs['volume']:.3f}")
>>>     if optimizer.converged():
>>>         break
__init__(problem: Problem, sub_tol=1e-07, sub_maxiter=100, change_tol=0.0001, fun_tol=1e-06, move=0.5, timer=False)[source]

Create an optimizer bound to a problem.

Parameters - problem: instance of Problem providing objective and gradients.

iter()[source]

Perform one MMA optimization iteration.

Updates design variables by solving a convex sub-problem approximation of the original problem. Uses moving asymptotes to control step size.

Returns:

If timer=True, returns iteration time in seconds

Return type:

float, optional

Notes

  • Solves convex sub-problem using mmasub()

  • Updates moving asymptotes (low, upp) for next iteration

  • Tracks design variable change and objective change

  • Handles NaN values by clipping to bounds

converged(*args, **kwargs)[source]

Check if optimizer has converged.

Returns:

True if convergence criteria are met: - Problem penalty continuation is complete (is_terminal() == True) - Design variable change <= change_tol - Objective function change <= fun_tol

Return type:

bool

Notes

Convergence requires both design change and objective change to be below tolerances. Also checks that penalty continuation (if used) has finished.

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘objective’: Current objective function value - ‘variable change’: L2 norm of design variable change - ‘function change’: Relative objective function change - Additional keys from problem.logs() (e.g., ‘iteration’, ‘residual’)

Return type:

dict

Notes

Used for monitoring optimization progress and convergence.

class pyFANTOM.Optimizers.CPU.OC.OC(problem: Problem, move: float = 0.2, change_tol=0.0001, fun_tol=1e-06, timer=False)[source]

Bases: Optimizer

Optimality Criteria (OC) method for topology optimization.

Classic heuristic optimizer specifically designed for topology optimization with independent volume constraints. Very fast but limited to specific problem formulations.

Parameters:

  • problem (Problem) – Optimization problem (must have independent constraints)

  • move (float, optional) – Move limit for design variables (default: 0.2)

  • change_tol (float, optional) – Convergence tolerance for design change (default: 1e-4)

  • fun_tol (float, optional) – Convergence tolerance for objective change (default: 1e-6)

  • timer (bool, optional) – Return iteration time (default: False)

ocP

OC update parameter computed from sensitivities

Type:

ndarray

iter()[source]

Perform one OC iteration

converged()[source]

Check convergence

Notes

  • Best for: Classic minimum compliance with volume constraint

  • Fastest: 2-3x faster than MMA per iteration

  • Limited scope: Only works with independent constraints

  • Heuristic: Not guaranteed to find KKT point

  • Typical use: Quick prototyping, benchmarking

  • Update formula: x_new = x * sqrt(-df/dg) projected onto move limits

  • Bisection: Finds Lagrange multiplier via bisection on volume constraint

Examples

>>> from pyFANTOM.CPU import OC
>>> optimizer = OC(problem=problem, move=0.2)
>>> for i in range(200):
>>>     optimizer.iter()
>>>     print(f"Iter {i}: C={optimizer.problem.f():.2e}")
>>>     if optimizer.converged():
>>>         break
__init__(problem: Problem, move: float = 0.2, change_tol=0.0001, fun_tol=1e-06, timer=False)[source]

Create an optimizer bound to a problem.

Parameters - problem: instance of Problem providing objective and gradients.

iter()[source]

Perform one OC optimization iteration.

Updates design variables using optimality criteria update formula with bisection to find Lagrange multiplier satisfying volume constraint.

Returns:

If timer=True, returns iteration time in seconds

Return type:

float, optional

Notes

  • Computes OC parameter: ocP = x * sqrt(-df/dg)

  • Uses bisection to find Lagrange multiplier lambda

  • Updates: x_new = clip(xL, xU, ocP/lambda)

  • Enforces move limits: xL = x - move, xU = x + move

converged(*args, **kwargs)[source]

Check if optimizer has converged.

Returns:

True if convergence criteria are met: - Problem penalty continuation is complete (is_terminal() == True) - Design variable change <= change_tol - Objective function change <= fun_tol

Return type:

bool

Notes

Convergence requires both design change and objective change to be below tolerances. Also checks that penalty continuation (if used) has finished.

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘objective’: Current objective function value - ‘variable change’: L2 norm of design variable change - ‘function change’: Relative objective function change - Additional keys from problem.logs() (e.g., ‘iteration’, ‘residual’)

Return type:

dict

Notes

Used for monitoring optimization progress and convergence.

class pyFANTOM.Optimizers.CPU.PGD.PGD(problem: Problem, change_tol=0.0001, fun_tol=1e-06, maxiter_N=50, tol_B=1e-08, tol_N=1e-06, C=1000000000000.0, fall_back_move=0.2, alpha_max=100.0, relaxation=1.0, warmup_iter=50, timer=False, conjugate_directions=True)[source]

Bases: Optimizer

Projected Gradient Descent optimizer with adaptive step sizes.

Fast first-order optimizer using BFGS-based step size adaptation and projection onto feasible set. Simpler and faster than MMA but less robust for complex constraints.

Parameters:

  • problem (Problem) – Optimization problem

  • change_tol (float, optional) – Design variable change tolerance (default: 1e-4)

  • fun_tol (float, optional) – Objective function change tolerance (default: 1e-6)

  • maxiter_N (int, optional) – Max iterations for Newton’s method in projection (default: 50)

  • tol_B (float, optional) – Bisection tolerance (default: 1e-8)

  • tol_N (float, optional) – Newton solver tolerance (default: 1e-6)

  • C (float, optional) – Large constant for penalty (default: 1e12)

  • fall_back_move (float, optional) – Fallback step size when BFGS fails (default: 0.2)

  • alpha_max (float, optional) – Maximum step size (default: 1e2)

  • relaxation (float, optional) – Step size relaxation parameter (default: 1.0)

  • warmup_iter (int, optional) – Iterations before enforcing strict feasibility (default: 50)

  • timer (bool, optional) – Return iteration time (default: False)

  • conjugate_directions (bool, optional) – Use conjugate gradient directions (default: True)

Notes

  • Best for: Simple bounds, single constraint problems

  • Faster than MMA: ~30-50% faster per iteration

  • Less robust: May struggle with multiple nonlinear constraints

  • BFGS step size: Adapts automatically from gradient history

  • Conjugate directions: Accelerates convergence

Examples

>>> from pyFANTOM.CPU import PGD
>>> optimizer = PGD(problem=problem, conjugate_directions=True)
>>> for i in range(100):
>>>     optimizer.iter()
>>>     if optimizer.converged():
>>>         break
__init__(problem: Problem, change_tol=0.0001, fun_tol=1e-06, maxiter_N=50, tol_B=1e-08, tol_N=1e-06, C=1000000000000.0, fall_back_move=0.2, alpha_max=100.0, relaxation=1.0, warmup_iter=50, timer=False, conjugate_directions=True)[source]

Create an optimizer bound to a problem.

Parameters - problem: instance of Problem providing objective and gradients.

alpha(desvars, df)[source]
project_to_feasible(desvars_new, desvars, dg)[source]
iter()[source]

Perform one PGD optimization iteration.

Updates design variables using projected gradient descent with adaptive step size (BFGS-based) and projection onto feasible set.

Returns:

If timer=True, returns iteration time in seconds

Return type:

float, optional

Notes

  • Computes step size using BFGS approximation (alpha)

  • Uses conjugate gradient directions if enabled

  • Projects update onto feasible set (satisfies constraints)

  • For independent constraints: uses bisection

  • For coupled constraints: uses Newton’s method with Wolfe line search

converged(*args, **kwargs)[source]

Check if optimizer has converged.

Returns:

True if convergence criteria are met: - Problem penalty continuation is complete (is_terminal() == True) - Design variable change <= change_tol - Objective function change <= fun_tol

Return type:

bool

Notes

Convergence requires both design change and objective change to be below tolerances. Also checks that penalty continuation (if used) has finished.

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘objective’: Current objective function value - ‘variable change’: L2 norm of design variable change - ‘function change’: Relative objective function change - Additional keys from problem.logs() (e.g., ‘iteration’, ‘residual’)

Return type:

dict

Notes

Used for monitoring optimization progress and convergence.

Wolfe conditions line search for Newton step

class pyFANTOM.Problem.CPU.MinimumCompliance.MinimumCompliance(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, E_mul: list[float] = [1.0], void: float = 1e-06, penalty: float = 3.0, volume_fraction: list[float] = [0.25], penalty_schedule: Callable[[float, int], float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5)[source]

Bases: Problem

Minimum compliance topology optimization problem.

Classic topology optimization: minimize structure compliance (maximize stiffness) subject to volume constraint. Implements SIMP (Solid Isotropic Material with Penalization) with optional density filtering and Heaviside projection for manufacturability.

Parameters:

  • FE (FiniteElement) – Finite element analysis engine with mesh, kernel, and solver

  • filter (StructuredFilter2D, StructuredFilter3D, or GeneralFilter) – Density filter for ensuring minimum feature sizes

  • E_mul (list of float, optional) – Young’s modulus multipliers for each material (default: [1.0] for single material)

  • void (float, optional) – Minimum density to avoid singularity (default: 1e-6)

  • penalty (float, optional) – SIMP penalization exponent (default: 3.0). Higher = more binary designs

  • volume_fraction (list of float, optional) – Volume fraction constraint for each material (default: [0.25])

  • penalty_schedule (callable, optional) – Function(p, iteration) for penalty continuation. If None, uses constant penalty

  • heavyside (bool, optional) – Apply Heaviside projection for sharper 0-1 designs (default: True)

  • beta (float or callable, optional) – Heaviside projection sharpness parameter. Can be a float (default: 2) or a callable function of iteration: beta(iteration) -> float. Enables beta continuation for gradual Heaviside sharpening during optimization.

  • eta (float, optional) – Heaviside projection threshold (default: 0.5)

is_single_material

True for single material, False for multi-material optimization

Type:

bool

n_material

Number of materials

Type:

int

iteration

Current iteration count

Type:

int

f()[source]

Compute compliance objective

nabla_f()[source]

Compute compliance sensitivities

g()[source]

Compute volume constraint(s)

nabla_g()[source]

Compute volume constraint gradients

get_desvars()[source]

Get current design variables (filtered densities)

set_desvars(rho)[source]

Set design variables and trigger FEA solve

visualize_solution(**kwargs)[source]

Plot optimized design

ill_conditioned()[source]

Check if FEA residual indicates poor conditioning

Notes

SIMP Penalization: - E_effective = E_min + rho^p * (E - E_min) - p=1: Linear, p=3: Standard, p>3: More binary - Use penalty_schedule for continuation: start low, increase gradually

Heaviside Projection: - Smooths 0-1 transition for manufacturing - beta controls sharpness (higher = sharper) - beta can be a float or callable(iteration) for continuation schedules - eta controls threshold location (0.5 = centered)

Multi-Material: - E_mul = [1.0, 0.5] creates two materials with different stiffnesses - volume_fraction = [0.3, 0.2] enforces separate volume constraints - Design variables shape: (n_materials * n_elements,)

Sensitivity Computation: - Uses adjoint method: dC/drho = -U^T @ dK/drho @ U - Filter adjoint applied via filter._rmatvec()

Examples

>>> from pyFANTOM.CPU import *
>>> # Setup FEA
>>> mesh = StructuredMesh2D(nx=128, ny=64, lx=2.0, ly=1.0)
>>> kernel = StructuredStiffnessKernel(mesh=mesh)
>>> solver = CHOLMOD(kernel=kernel)
>>> FE = FiniteElement(mesh=mesh, kernel=kernel, solver=solver)
>>>
>>> # Apply BCs and loads
>>> FE.add_dirichlet_boundary_condition(node_ids=fixed_nodes, rhs=0)
>>> FE.add_point_forces(node_ids=load_nodes, forces=np.array([[0, -1.0]]))
>>>
>>> # Define optimization problem
>>> filter = StructuredFilter2D(mesh=mesh, r_min=1.5)
>>> problem = MinimumCompliance(FE=FE, filter=filter, penalty=3.0, volume_fraction=[0.3])
>>>
>>> # Multi-material example
>>> problem_multi = MinimumCompliance(
>>>     FE=FE, filter=filter,
>>>     E_mul=[1.0, 0.5],  # Two materials
>>>     volume_fraction=[0.25, 0.15],  # Separate constraints
>>>     penalty=3.0
>>> )
__init__(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, E_mul: list[float] = [1.0], void: float = 1e-06, penalty: float = 3.0, volume_fraction: list[float] = [0.25], penalty_schedule: Callable[[float, int], float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5)[source]

Initialize problem state.

Subclasses may accept finite-element handlers, filters, and other configuration parameters.

N()[source]

Return the number of design variables.

Returns:

Total number of design variables (n_elements * n_materials)

Return type:

int

m()[source]

Return the number of constraints.

Returns:

Number of volume constraints (1 for single material, n_materials for multi-material)

Return type:

int

is_independent()[source]

Check if constraints are independent (required for OC optimizer).

Returns:

True if constraints are independent (always True for MinimumCompliance)

Return type:

bool

constraint_map()[source]

Return mapping of constraints to design variables.

Returns:

  • Single material: 1 (scalar)

  • Multi-material: array of shape (n_materials, n_vars) with 1s indicating which design variables belong to each material constraint

Return type:

int or ndarray

Notes

Used by optimizers to identify which design variables affect each constraint. For multi-material, constraint i affects variables [i*n_elements:(i+1)*n_elements].

bounds()[source]

Return bounds for design variables.

Returns:

(lower_bound, upper_bound) where both are 0.0 and 1.0 respectively

Return type:

tuple

Notes

Design variables represent normalized densities in [0, 1].

visualize_problem(**kwargs)[source]

Visualize problem setup (mesh, BCs, loads).

Parameters:

**kwargs – Arguments passed to FE.visualize_problem()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

visualize_solution(**kwargs)[source]

Visualize optimized design (density distribution).

Parameters:

**kwargs – Arguments passed to FE.visualize_density()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Shows current design variables as density field. For multi-material problems, displays combined material distribution.

init_desvars()[source]

Initialize design variables to uniform density at volume fraction.

Sets all design variables to the volume fraction constraint value and performs initial FEA solve to compute objective and constraints.

Notes

  • Single material: all variables set to volume_fraction

  • Multi-material: all variables set to min(volume_fraction) for each material

  • Resets iteration counter to 0

  • Triggers _compute() to evaluate objective and constraints

set_desvars(desvars: ndarray)[source]

Set design variables and recompute objective/constraints.

Parameters:

desvars (ndarray) – Design variables, shape (n_vars,). Values should be in [0, 1]

Raises:

ValueError – If desvars shape doesn’t match num_vars

Notes

  • Triggers FEA solve and sensitivity computation

  • Increments iteration counter

  • Updates cached values: _f, _g, _nabla_f, _nabla_g

get_desvars()[source]

Get current design variables.

Returns:

Current design variables, shape (n_vars,)

Return type:

ndarray

Notes

Returns raw (unfiltered) design variables. For filtered densities, access via problem.filter.dot(problem.get_desvars()).

penalize(rho: ndarray)[source]
penalize_grad(rho: ndarray)[source]
f(rho: ndarray | None = None)[source]

Compute objective function value (compliance).

Parameters:

rho (ndarray, optional) – Design variables for linearization. If None, returns cached value

Returns:

Compliance value (F^T @ U). Lower is better (stiffer structure).

Return type:

float

Notes

  • If rho provided: returns linearized approximation f(x) + df/dx @ (rho - x)

  • If rho is None: returns cached value from last set_desvars() call

  • Compliance = F^T @ U = U^T @ K @ U (strain energy)

nabla_f(rho: ndarray | None = None)[source]

Compute objective function gradient (compliance sensitivities).

Parameters:

rho (ndarray, optional) – Unused (for interface compatibility)

Returns:

Gradient of compliance w.r.t. design variables, shape (n_vars,)

Return type:

ndarray

Notes

  • Uses adjoint method: dC/drho = -U^T @ dK/drho @ U

  • Includes filter adjoint: sens_raw = H^T @ sens_filtered

  • Negative gradient means increasing density reduces compliance (good)

g(rho=None)[source]

Compute constraint values (volume fraction violations).

Parameters:

rho (ndarray, optional) – Design variables. If None, uses current desvars

Returns:

Constraint violations, shape (n_constraints,). Negative = satisfied. g[i] = (volume_fraction[i] - actual_volume_fraction[i])

Return type:

ndarray

Notes

  • Constraint satisfied when g <= 0

  • For single material: returns scalar

  • For multi-material: returns array with one constraint per material

nabla_g()[source]

Compute constraint gradients (volume fraction sensitivities).

Returns:

Gradient of constraints w.r.t. design variables. - Single material: shape (n_vars,) - Multi-material: shape (n_materials, n_vars)

Return type:

ndarray

Notes

  • Each row is d(volume_fraction[i])/drho

  • For uniform elements, gradient is constant: 1/volume_total per element

  • Used by optimizers to enforce volume constraints

ill_conditioned()[source]

Check if FEA system is ill-conditioned.

Returns:

True if residual >= 1e-2 (indicates poor solver convergence)

Return type:

bool

Notes

  • Residual > 1e-2 suggests numerical issues (check void parameter, penalty)

  • May indicate near-singular stiffness matrix (too many void elements)

  • Consider increasing void parameter or reducing penalty

is_terminal()[source]

Check if penalty continuation has reached final value.

Returns:

True if penalty schedule is complete or not used

Return type:

bool

Notes

  • Used by optimizers to determine if continuation is finished

  • If penalty_schedule is None, always returns True

  • If penalty_schedule exists, checks if current iteration has reached final penalty

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘iteration’: Current iteration number - ‘residual’: FEA solver residual (||K@U - F|| / ||F||)

Return type:

dict

Notes

Used by optimizers to track convergence and diagnose issues.

FEA(thresshold: bool = True)[source]
visualize_field(field, ax=None, thresshold=True, **kwargs)[source]
class pyFANTOM.Problem.CPU.MinimumComplianceNL.MinimumComplianceNL(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, E_mul: list[float] = [1.0], void: float = 1e-06, penalty: float = 3.0, volume_fraction: list[float] = [0.25], penalty_schedule: list[float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5, independant: bool = True)[source]

Bases: Problem

Minimum compliance topology optimization for nonlinear elasticity.

Nonlinear version of MinimumCompliance using geometrically nonlinear finite element analysis. Handles large deformations and material nonlinearity via Newton-Raphson iterations. Requires NLFiniteElement and NLElasticity physics.

Parameters:

  • FE (NLFiniteElement) – Nonlinear finite element analysis engine with mesh, kernel, and solver

  • filter (StructuredFilter2D, StructuredFilter3D, or GeneralFilter) – Density filter for ensuring minimum feature sizes

  • E_mul (list of float, optional) – Young’s modulus multipliers for each material (default: [1.0] for single material)

  • void (float, optional) – Minimum density to avoid singularity (default: 1e-6)

  • penalty (float, optional) – SIMP penalization exponent (default: 3.0). Higher = more binary designs

  • volume_fraction (list of float, optional) – Volume fraction constraint for each material (default: [0.25])

  • penalty_schedule (callable, optional) – Function(p, iteration) for penalty continuation. If None, uses constant penalty

  • heavyside (bool, optional) – Apply Heaviside projection for sharper 0-1 designs (default: True)

  • beta (float or callable, optional) – Heaviside projection sharpness parameter. Can be a float (default: 2) or a callable function of iteration: beta(iteration) -> float. Enables beta continuation for gradual Heaviside sharpening during optimization.

  • eta (float, optional) – Heaviside projection threshold (default: 0.5)

  • independant (bool, optional) – Whether constraints are independent (default: True)

is_single_material

True for single material, False for multi-material optimization

Type:

bool

n_material

Number of materials

Type:

int

iteration

Current iteration count

Type:

int

f()[source]

Compute compliance objective

nabla_f()[source]

Compute compliance sensitivities (nonlinear adjoint)

g()[source]

Compute volume constraint(s)

nabla_g()[source]

Compute volume constraint gradients

get_desvars()[source]

Get current design variables (filtered densities)

set_desvars(rho)[source]

Set design variables and trigger nonlinear FEA solve

visualize_solution(**kwargs)[source]

Plot optimized design

Notes

Nonlinear Analysis: - Uses Newton-Raphson method for equilibrium iterations - Tangent stiffness matrix updated each iteration - Load stepping for convergence (4 load steps by default) - More computationally expensive than linear MinimumCompliance

Sensitivity Analysis: - Uses nonlinear adjoint method - Requires solving additional adjoint system - Sensitivities account for geometric nonlinearity

Use Cases: - Large deformation problems - Buckling-sensitive structures - Compliant mechanisms - Soft robotics applications

Examples

>>> from pyFANTOM.CPU import *
>>> from pyFANTOM import NLElasticity
>>> # Setup nonlinear FEA
>>> physics = NLElasticity(E=1.0, nu=0.3)
>>> mesh = StructuredMesh2D(nx=64, ny=32, lx=2.0, ly=1.0, physics=physics)
>>> kernel = NLUniformStiffnessKernel(mesh=mesh)
>>> solver = CG(kernel=kernel)
>>> FE = NLFiniteElement(mesh=mesh, kernel=kernel, solver=solver, physics=physics)
>>>
>>> # Apply BCs and loads
>>> FE.add_dirichlet_boundary_condition(node_ids=fixed_nodes, rhs=0)
>>> FE.add_point_forces(node_ids=load_nodes, forces=np.array([[0, -1.0]]))
>>>
>>> # Define optimization problem
>>> filter = StructuredFilter2D(mesh=mesh, r_min=1.5)
>>> problem = MinimumComplianceNL(FE=FE, filter=filter, penalty=3.0, volume_fraction=[0.3])
__init__(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, E_mul: list[float] = [1.0], void: float = 1e-06, penalty: float = 3.0, volume_fraction: list[float] = [0.25], penalty_schedule: list[float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5, independant: bool = True)[source]

Initialize problem state.

Subclasses may accept finite-element handlers, filters, and other configuration parameters.

N()[source]

Return the number of design variables.

Returns:

Total number of design variables (n_elements * n_materials)

Return type:

int

m()[source]

Return the number of constraints.

Returns:

Number of volume constraints (1 for single material, n_materials for multi-material)

Return type:

int

is_independent()[source]

Check if constraints are independent (required for OC optimizer).

Returns:

True if constraints are independent (controlled by independant parameter)

Return type:

bool

constraint_map()[source]

Return mapping of constraints to design variables.

Returns:

  • Single material: 1 (scalar)

  • Multi-material: array of shape (n_materials, n_vars) with 1s indicating which design variables belong to each material constraint

Return type:

int or ndarray

Notes

Used by optimizers to identify which design variables affect each constraint. For multi-material, constraint i affects variables [i*n_elements:(i+1)*n_elements].

bounds()[source]

Return bounds for design variables.

Returns:

(lower_bound, upper_bound) where both are 0.0 and 1.0 respectively

Return type:

tuple

Notes

Design variables represent normalized densities in [0, 1].

visualize_problem(**kwargs)[source]

Visualize problem setup (mesh, BCs, loads).

Parameters:

**kwargs – Arguments passed to FE.visualize_problem()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

visualize_solution(**kwargs)[source]

Visualize optimized design (density distribution).

Parameters:

**kwargs – Arguments passed to FE.visualize_density()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Shows current design variables as density field. For multi-material problems, displays combined material distribution.

init_desvars()[source]

Initialize design variables to uniform density at volume fraction.

Sets all design variables to the volume fraction constraint value and performs initial nonlinear FEA solve to compute objective and constraints.

Notes

  • Single material: all variables set to volume_fraction

  • Multi-material: all variables set to volume_fraction for each material

  • Resets iteration counter to 0

  • Triggers _compute() to evaluate objective and constraints via nonlinear solve

set_desvars(desvars: ndarray)[source]

Set design variables and recompute objective/constraints.

Parameters:

desvars (ndarray) – Design variables, shape (n_vars,). Values should be in [0, 1]

Raises:

ValueError – If desvars shape doesn’t match num_vars

Notes

  • Triggers nonlinear FEA solve (Newton-Raphson) and sensitivity computation

  • Increments iteration counter

  • Updates cached values: _f, _g, _nabla_f, _nabla_g

  • More expensive than linear MinimumCompliance due to nonlinear solve

get_desvars()[source]

Get current design variables.

Returns:

Current design variables, shape (n_vars,)

Return type:

ndarray

Notes

Returns raw (unfiltered) design variables. For filtered densities, access via problem.filter.dot(problem.get_desvars()).

penalize(rho: ndarray)[source]

Apply SIMP penalization to densities.

Parameters:

rho (ndarray) – Filtered design variables, shape (n_elements,)

Returns:

Penalized Young’s modulus values, shape (n_elements,)

Return type:

ndarray

Notes

Applies rho^penalty and scales by E_mul. Clips to void to avoid singularity. For nonlinear problems, uses simplified penalization (no Heaviside).

penalize_grad(rho: ndarray)[source]

Compute gradient of SIMP penalization.

Parameters:

rho (ndarray) – Filtered design variables, shape (n_elements,)

Returns:

Gradient of penalization w.r.t. rho, shape (n_elements,)

Return type:

ndarray

Notes

Computes d(rho^penalty * E_mul)/drho = penalty * rho^(penalty-1) * E_mul. Used in sensitivity analysis for nonlinear problems.

f(rho: ndarray | None = None)[source]

Compute objective function value (compliance).

Parameters:

rho (ndarray, optional) – Design variables for linearization. If None, returns cached value

Returns:

Compliance value (F^T @ U). Lower is better (stiffer structure).

Return type:

float

Notes

  • If rho provided: returns linearized approximation f(x) + df/dx @ (rho - x)

  • If rho is None: returns cached value from last set_desvars() call

  • Compliance computed from nonlinear equilibrium solution

nabla_f(rho: ndarray | None = None)[source]

Compute objective function gradient (compliance sensitivities).

Parameters:

rho (ndarray, optional) – Unused (for interface compatibility)

Returns:

Gradient of compliance w.r.t. design variables, shape (n_vars,)

Return type:

ndarray

Notes

  • Uses nonlinear adjoint method: dC/drho = lambda^T @ dR/drho

  • Includes filter adjoint: sens_raw = H^T @ sens_filtered

  • Accounts for geometric nonlinearity in sensitivity computation

g(rho=None)[source]

Compute constraint values (volume fraction violations).

Parameters:

rho (ndarray, optional) – Design variables. If None, uses current desvars

Returns:

Constraint violations, shape (n_constraints,). Negative = satisfied. g[i] = (volume_fraction[i] - actual_volume_fraction[i])

Return type:

ndarray

Notes

  • Constraint satisfied when g <= 0

  • For single material: returns scalar

  • For multi-material: returns array with one constraint per material

nabla_g()[source]

Compute constraint gradients (volume fraction sensitivities).

Returns:

Gradient of constraints w.r.t. design variables. - Single material: shape (n_vars,) - Multi-material: shape (n_materials, n_vars)

Return type:

ndarray

Notes

  • Each row is d(volume_fraction[i])/drho

  • For uniform elements, gradient is constant: 1/volume_total per element

  • Used by optimizers to enforce volume constraints

ill_conditioned()[source]

Check if FEA system is ill-conditioned.

Returns:

True if residual >= 1e-2 (indicates poor solver convergence)

Return type:

bool

Notes

  • Residual > 1e-2 suggests numerical issues (check void parameter, penalty)

  • May indicate near-singular tangent stiffness matrix

  • Consider increasing void parameter or reducing penalty

  • For nonlinear problems, may also indicate convergence issues in Newton-Raphson

is_terminal()[source]

Check if penalty continuation has reached final value.

Returns:

True if penalty schedule is complete or not used

Return type:

bool

Notes

  • Used by optimizers to determine if continuation is finished

  • If penalty_schedule is None, always returns True

  • If penalty_schedule exists, checks if current iteration has reached final penalty

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘iteration’: Current iteration number - ‘residual’: FEA solver residual (||K@U - F|| / ||F||)

Return type:

dict

Notes

Used by optimizers to track convergence and diagnose issues. Residual includes both linear solver and Newton-Raphson convergence.

FEA(thresshold: bool = True)[source]
visualize_field(field, ax=None, thresshold=True, **kwargs)[source]
class pyFANTOM.Problem.CPU.ComplianceConstrainedMinimumVolume.ComplianceConstrainedMinimumVolume(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, E: float = 1.0, void: float = 1e-06, penalty: float = 3.0, compliance_limit: float = 0.25, penalty_schedule: Callable[[float, int], float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5)[source]

Bases: Problem

Minimum volume topology optimization with compliance constraint.

Inverse formulation of MinimumCompliance: minimize material usage while maintaining structure stiffness above threshold. Useful for lightweight design with performance requirements.

Parameters:

  • FE (FiniteElement) – Finite element analysis engine

  • filter (StructuredFilter2D, StructuredFilter3D, or GeneralFilter) – Density filter

  • E (float, optional) – Young’s modulus (default: 1.0)

  • void (float, optional) – Minimum density to avoid singularity (default: 1e-6)

  • penalty (float, optional) – SIMP penalization exponent (default: 3.0)

  • compliance_limit (float, optional) – Maximum allowable compliance (normalized by initial compliance) (default: 0.25)

  • penalty_schedule (callable, optional) – Penalty continuation function(p, iteration)

  • heavyside (bool, optional) – Apply Heaviside projection (default: True)

  • beta (float or callable, optional) – Heaviside projection sharpness parameter. Can be a float (default: 2) or a callable function of iteration: beta(iteration) -> float. Enables beta continuation for gradual Heaviside sharpening during optimization.

  • eta (float, optional) – Heaviside threshold (default: 0.5)

Notes

  • Objective: Minimize volume

  • Constraint: Compliance ≤ compliance_limit * initial_compliance

  • Use case: Lightweight design, material cost minimization

  • Comparison: More challenging than MinimumCompliance (compliance constraint is nonlinear)

  • Typical compliance_limit: 0.2-0.5 of initial design

Examples

>>> problem = ComplianceConstrainedMinimumVolume(
>>>     FE=FE, filter=filter,
>>>     compliance_limit=0.3,  # 30% of initial compliance
>>>     penalty=3.0
>>> )
__init__(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, E: float = 1.0, void: float = 1e-06, penalty: float = 3.0, compliance_limit: float = 0.25, penalty_schedule: Callable[[float, int], float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5)[source]

Initialize problem state.

Subclasses may accept finite-element handlers, filters, and other configuration parameters.

N()[source]

Return the number of design variables.

Returns:

Total number of design variables (n_elements)

Return type:

int

m()[source]

Return the number of constraints.

Returns:

Number of constraints (1: compliance constraint)

Return type:

int

is_independent()[source]

Check if constraints are independent (required for OC optimizer).

Returns:

True if constraints are independent (always True for ComplianceConstrainedMinimumVolume)

Return type:

bool

constraint_map()[source]

Return mapping of constraints to design variables.

Returns:

Scalar value 1 indicating all design variables affect the compliance constraint

Return type:

int

Notes

Used by optimizers to identify which design variables affect the constraint. For this problem, all variables affect compliance.

bounds()[source]

Return bounds for design variables.

Returns:

(lower_bound, upper_bound) where both are 0.0 and 1.0 respectively

Return type:

tuple

Notes

Design variables represent normalized densities in [0, 1].

visualize_problem(**kwargs)[source]

Visualize problem setup (mesh, BCs, loads).

Parameters:

**kwargs – Arguments passed to FE.visualize_problem()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

visualize_solution(**kwargs)[source]

Visualize optimized design (density distribution).

Parameters:

**kwargs – Arguments passed to FE.visualize_density()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Shows current design variables as density field.

init_desvars()[source]

Initialize design variables to full density (volume = 1.0).

Sets all design variables to 1.0 (full material) and performs initial FEA solve to compute objective and constraints. This ensures initial design satisfies compliance constraint.

Notes

  • All variables set to 1.0 (full material)

  • Resets iteration counter to 0

  • Triggers _compute() to evaluate objective and constraints

set_desvars(desvars: ndarray)[source]

Set design variables and recompute objective/constraints.

Parameters:

desvars (ndarray) – Design variables, shape (n_vars,). Values should be in [0, 1]

Raises:

ValueError – If desvars shape doesn’t match num_vars

Notes

  • Triggers FEA solve and sensitivity computation

  • Increments iteration counter

  • Updates cached values: _f, _g, _nabla_f, _nabla_g

get_desvars()[source]

Get current design variables.

Returns:

Current design variables, shape (n_vars,)

Return type:

ndarray

Notes

Returns raw (unfiltered) design variables. For filtered densities, access via problem.filter.dot(problem.get_desvars()).

penalize(rho: ndarray)[source]

Apply SIMP penalization and Heaviside projection.

Parameters:

rho (ndarray) – Filtered densities, shape (n_elements,)

Returns:

Penalized Young’s moduli, shape (n_elements,)

Return type:

ndarray

Notes

  • Applies Heaviside projection if enabled

  • Applies SIMP: E = E0 * rho^penalty

  • Clips to void value to avoid singularity

  • Uses penalty_schedule if provided

  • Uses beta_schedule if beta is callable

penalize_grad(rho: ndarray)[source]

Compute gradient of penalization function.

Parameters:

rho (ndarray) – Filtered densities, shape (n_elements,)

Returns:

Gradient dE/drho, shape (n_elements,)

Return type:

ndarray

Notes

  • Derivative of SIMP + Heaviside projection

  • Used in chain rule for sensitivity analysis

  • Accounts for penalty_schedule if provided

  • Accounts for beta_schedule if beta is callable

f(rho: ndarray | None = None)[source]

Compute objective function value (volume).

Parameters:

rho (ndarray, optional) – Design variables for linearization. If None, returns cached value

Returns:

Volume fraction (normalized by total domain volume). Lower is better.

Return type:

float

Notes

  • If rho provided: returns linearized approximation f(x) + df/dx @ (rho - x)

  • If rho is None: returns cached value from last set_desvars() call

  • Objective is to minimize material usage while satisfying compliance constraint

nabla_f(rho: ndarray | None = None)[source]

Compute objective function gradient (volume sensitivities).

Parameters:

rho (ndarray, optional) – Unused (for interface compatibility)

Returns:

Gradient of volume w.r.t. design variables, shape (n_vars,)

Return type:

ndarray

Notes

  • Constant gradient: 1/volume_total per element

  • All elements contribute equally to volume

  • Used by optimizers to minimize material usage

g(rho=None)[source]

Compute constraint values (compliance constraint violations).

Parameters:

rho (ndarray, optional) – Design variables. If None, uses current desvars

Returns:

Normalized compliance constraint violation. Negative = satisfied. g = (compliance - compliance_limit) / compliance_limit

Return type:

float

Notes

  • Constraint satisfied when g <= 0

  • Normalized by compliance_limit for better scaling

  • If rho provided, uses linearized approximation

nabla_g()[source]

Compute constraint gradients (compliance sensitivities).

Returns:

Gradient of compliance constraint w.r.t. design variables, shape (n_vars,)

Return type:

ndarray

Notes

  • Normalized by compliance_limit for better scaling

  • Uses adjoint method: dC/drho = -U^T @ dK/drho @ U

  • Includes filter adjoint: sens_raw = H^T @ sens_filtered

  • Used by optimizers to enforce compliance constraint

ill_conditioned()[source]

Check if FEA system is ill-conditioned.

Returns:

True if residual >= 1e-2 (indicates poor solver convergence)

Return type:

bool

Notes

  • Residual > 1e-2 suggests numerical issues (check void parameter, penalty)

  • May indicate near-singular stiffness matrix (too many void elements)

  • Consider increasing void parameter or reducing penalty

is_terminal()[source]

Check if penalty continuation has reached final value.

Returns:

True if penalty schedule is complete or not used

Return type:

bool

Notes

  • Used by optimizers to determine if continuation is finished

  • If penalty_schedule is None, always returns True

  • If penalty_schedule exists, checks if current iteration has reached final penalty

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘iteration’: Current iteration number - ‘residual’: FEA solver residual (||K@U - F|| / ||F||)

Return type:

dict

Notes

Used by optimizers to track convergence and diagnose issues.

class pyFANTOM.Problem.CPU.WeightDistributionMinimumCompliance.WeightDistributionMinimumCompliance(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, target_centroid: List[float], maximum_distance: float, E_mul: list[float] = [1.0], void: float = 1e-06, penalty: float = 3.0, volume_fraction: list[float] = [0.25], penalty_schedule: Callable[[float, int], float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5)[source]

Bases: Problem

Minimum compliance with weight distribution constraint.

Extends MinimumCompliance with an additional constraint on the centroid location of the material distribution. Useful for balancing structures or controlling center of mass for dynamic applications.

Parameters:

  • FE (FiniteElement) – Finite element analysis engine with mesh, kernel, and solver

  • filter (StructuredFilter2D, StructuredFilter3D, or GeneralFilter) – Density filter for ensuring minimum feature sizes

  • target_centroid (List[float]) – Target centroid location in physical coordinates, shape (2,) for 2D or (3,) for 3D

  • maximum_distance (float) – Maximum allowed distance from target centroid (constraint: ||centroid - target|| <= maximum_distance)

  • E_mul (list of float, optional) – Young’s modulus multipliers for each material (default: [1.0] for single material)

  • void (float, optional) – Minimum density to avoid singularity (default: 1e-6)

  • penalty (float, optional) – SIMP penalization exponent (default: 3.0). Higher = more binary designs

  • volume_fraction (list of float, optional) – Volume fraction constraint for each material (default: [0.25])

  • penalty_schedule (callable, optional) – Function(p, iteration) for penalty continuation. If None, uses constant penalty

  • heavyside (bool, optional) – Apply Heaviside projection for sharper 0-1 designs (default: True)

  • beta (float or callable, optional) – Heaviside projection sharpness parameter. Can be a float (default: 2) or a callable function of iteration: beta(iteration) -> float. Enables beta continuation for gradual Heaviside sharpening during optimization.

  • eta (float, optional) – Heaviside projection threshold (default: 0.5)

is_single_material

True for single material, False for multi-material optimization

Type:

bool

n_material

Number of materials

Type:

int

iteration

Current iteration count

Type:

int

target_centroid

Target centroid location

Type:

ndarray

maximum_distance

Maximum allowed distance from target

Type:

float

f()[source]

Compute compliance objective

nabla_f()[source]

Compute compliance sensitivities

g()[source]

Compute constraint values (volume + centroid distance)

nabla_g()[source]

Compute constraint gradients

get_desvars()[source]

Get current design variables

set_desvars(rho)[source]

Set design variables and trigger FEA solve

visualize_solution(**kwargs)[source]

Plot optimized design

Notes

Constraints: - Volume constraint: Same as MinimumCompliance - Centroid constraint: ||centroid(rho) - target_centroid|| <= maximum_distance - Constraints are coupled (not independent), so OC optimizer cannot be used

Centroid Computation: - Centroid = sum(rho[i] * volume[i] * centroid[i]) / sum(rho[i] * volume[i]) - Weighted by element density and volume

Use Cases: - Balancing structures for dynamic loads - Controlling center of mass - Symmetric designs with off-center loading - Multi-material designs with weight distribution requirements

Examples

>>> from pyFANTOM.CPU import *
>>> # Setup FEA
>>> mesh = StructuredMesh2D(nx=128, ny=64, lx=2.0, ly=1.0)
>>> kernel = StructuredStiffnessKernel(mesh=mesh)
>>> solver = CHOLMOD(kernel=kernel)
>>> FE = FiniteElement(mesh=mesh, kernel=kernel, solver=solver)
>>>
>>> # Apply BCs and loads
>>> FE.add_dirichlet_boundary_condition(node_ids=fixed_nodes, rhs=0)
>>> FE.add_point_forces(node_ids=load_nodes, forces=np.array([[0, -1.0]]))
>>>
>>> # Define optimization problem with centroid constraint
>>> filter = StructuredFilter2D(mesh=mesh, r_min=1.5)
>>> problem = WeightDistributionMinimumCompliance(
>>>     FE=FE, filter=filter,
>>>     target_centroid=[1.0, 0.5],  # Center of 2D domain
>>>     maximum_distance=0.1,  # 10% of domain size
>>>     penalty=3.0, volume_fraction=[0.3]
>>> )
__init__(FE: FiniteElement, filter: StructuredFilter2D | StructuredFilter3D | GeneralFilter, target_centroid: List[float], maximum_distance: float, E_mul: list[float] = [1.0], void: float = 1e-06, penalty: float = 3.0, volume_fraction: list[float] = [0.25], penalty_schedule: Callable[[float, int], float] | None = None, heavyside: bool = True, beta: float | Callable[[int], float] = 2, eta: float = 0.5)[source]

Initialize problem state.

Subclasses may accept finite-element handlers, filters, and other configuration parameters.

N()[source]

Return the number of design variables.

Returns:

Total number of design variables (n_elements * n_materials)

Return type:

int

m()[source]

Return the number of constraints.

Returns:

Number of constraints: - Single material: 2 (volume + centroid distance) - Multi-material: n_materials + 1 (volume per material + centroid distance)

Return type:

int

is_independent()[source]

Check if constraints are independent (required for OC optimizer).

Returns:

Always False. Centroid constraint couples all design variables, making constraints non-independent. Use MMA or PGD optimizer instead of OC.

Return type:

bool

constraint_map()[source]

Return mapping of constraints to design variables.

Returns:

  • Single material: 1 (scalar, all variables affect both constraints)

  • Multi-material: array of shape (n_materials + 1, n_vars) where: - Rows 0 to n_materials-1: volume constraints per material - Last row: centroid constraint (all variables = 1)

Return type:

int or ndarray

Notes

Used by optimizers to identify which design variables affect each constraint. The centroid constraint depends on all design variables, making constraints coupled.

bounds()[source]

Return bounds for design variables.

Returns:

(lower_bound, upper_bound) where both are 0.0 and 1.0 respectively

Return type:

tuple

Notes

Design variables represent normalized densities in [0, 1].

visualize_problem(**kwargs)[source]

Visualize problem setup (mesh, BCs, loads).

Parameters:

**kwargs – Arguments passed to FE.visualize_problem()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

visualize_solution(**kwargs)[source]

Visualize optimized design (density distribution).

Parameters:

**kwargs – Arguments passed to FE.visualize_density()

Returns:

Plot object

Return type:

matplotlib.axes.Axes or k3d.Plot

Notes

Shows current design variables as density field. For multi-material problems, displays combined material distribution.

init_desvars()[source]

Initialize design variables with random smoothed distribution.

Uses random initialization followed by triple filtering to create a smooth initial guess. This helps avoid getting stuck in poor local minima that may violate the centroid constraint.

Notes

  • Random initialization: uniform distribution in [min(0, 2*vf-1), 1.0]

  • Triple filtering: applies filter 3 times for smoothness

  • Normalized to match volume fraction constraint

  • Resets iteration counter to 0

  • Triggers _compute() to evaluate objective and constraints

set_desvars(desvars: ndarray)[source]

Set design variables and recompute objective/constraints.

Parameters:

desvars (ndarray) – Design variables, shape (n_vars,). Values should be in [0, 1]

Raises:

ValueError – If desvars shape doesn’t match num_vars

Notes

  • Triggers FEA solve and sensitivity computation

  • Increments iteration counter

  • Updates cached values: _f, _g, _nabla_f, _nabla_g

  • Computes centroid and centroid constraint gradient

get_desvars()[source]

Get current design variables.

Returns:

Current design variables, shape (n_vars,)

Return type:

ndarray

Notes

Returns raw (unfiltered) design variables. For filtered densities, access via problem.filter.dot(problem.get_desvars()).

penalize(rho: ndarray)[source]

Apply SIMP penalization and Heaviside projection.

Parameters:

rho (ndarray) – Filtered densities, shape (n_elements,) or (n_elements, n_materials)

Returns:

Penalized Young’s moduli, shape (n_elements,)

Return type:

ndarray

Notes

  • Applies Heaviside projection if enabled

  • Applies SIMP: E = E0 * rho^penalty

  • Clips to void value to avoid singularity

  • Uses penalty_schedule if provided

  • Uses beta_schedule if beta is callable

  • For multi-material: uses material interpolation scheme

penalize_grad(rho: ndarray)[source]

Compute gradient of penalization function.

Parameters:

rho (ndarray) – Filtered densities, shape (n_elements,) or (n_elements, n_materials)

Returns:

Gradient dE/drho, shape (n_elements,) or (n_elements, n_materials)

Return type:

ndarray

Notes

  • Derivative of SIMP + Heaviside projection

  • Used in chain rule for sensitivity analysis

  • Accounts for penalty_schedule if provided

  • Accounts for beta_schedule if beta is callable

  • For multi-material: returns gradient per material

f(rho: ndarray | None = None)[source]

Compute objective function value (compliance).

Parameters:

rho (ndarray, optional) – Design variables for linearization. If None, returns cached value

Returns:

Compliance value (F^T @ U). Lower is better (stiffer structure).

Return type:

float

Notes

  • If rho provided: returns linearized approximation f(x) + df/dx @ (rho - x)

  • If rho is None: returns cached value from last set_desvars() call

  • Compliance = F^T @ U = U^T @ K @ U (strain energy)

nabla_f(rho: ndarray | None = None)[source]

Compute objective function gradient (compliance sensitivities).

Parameters:

rho (ndarray, optional) – Unused (for interface compatibility)

Returns:

Gradient of compliance w.r.t. design variables, shape (n_vars,)

Return type:

ndarray

Notes

  • Uses adjoint method: dC/drho = -U^T @ dK/drho @ U

  • Includes filter adjoint: sens_raw = H^T @ sens_filtered

  • Negative gradient means increasing density reduces compliance (good)

g(rho=None)[source]

Compute constraint values (volume + centroid distance violations).

Parameters:

rho (ndarray, optional) – Design variables. If None, uses current desvars

Returns:

Constraint violations, shape (n_constraints,). Negative = satisfied. - First n_materials constraints: volume fraction violations - Last constraint: centroid distance violation (distance^2 - maximum_distance^2)

Return type:

ndarray

Notes

  • Constraint satisfied when g <= 0

  • For single material: returns [volume_violation, centroid_violation]

  • For multi-material: returns [vol_viol_1, …, vol_viol_n, centroid_violation]

  • If rho provided, uses linearized approximation

nabla_g()[source]

Compute constraint gradients (volume + centroid sensitivities).

Returns:

Gradient of constraints w.r.t. design variables. - Single material: shape (2, n_vars) - [volume_grad, centroid_grad] - Multi-material: shape (n_materials + 1, n_vars) - [vol_grad_1, …, vol_grad_n, centroid_grad]

Return type:

ndarray

Notes

  • First rows: d(volume_fraction[i])/drho (constant for uniform elements)

  • Last row: d(centroid_distance^2)/drho (depends on element positions)

  • Centroid gradient accounts for weighted average of element centroids

  • Used by optimizers to enforce constraints

ill_conditioned()[source]

Check if FEA system is ill-conditioned.

Returns:

True if residual >= 1e-2 (indicates poor solver convergence)

Return type:

bool

Notes

  • Residual > 1e-2 suggests numerical issues (check void parameter, penalty)

  • May indicate near-singular stiffness matrix (too many void elements)

  • Consider increasing void parameter or reducing penalty

is_terminal()[source]

Check if penalty continuation has reached final value.

Returns:

True if penalty schedule is complete or not used

Return type:

bool

Notes

  • Used by optimizers to determine if continuation is finished

  • If penalty_schedule is None, always returns True

  • If penalty_schedule exists, checks if current iteration has reached final penalty

logs()[source]

Return diagnostic information for current iteration.

Returns:

Dictionary with keys: - ‘iteration’: Current iteration number - ‘residual’: FEA solver residual (||K@U - F|| / ||F||)

Return type:

dict

Notes

Used by optimizers to track convergence and diagnose issues.

FEA(thresshold: bool = True)[source]
visualize_field(field, ax=None, thresshold=True, **kwargs)[source]