sktopt.core

class sktopt.core.DensityMethod(cfg: DensityMethodConfig, tsk: FEMDomain)

Bases: DensityMethodBase

Core driver for sensitivity-based density topology optimization.

This class implements the common workflow shared by multiple optimization strategies (e.g., OC, MOC). It handles

• problem setup (scaling, restart, output directories),

• finite element analysis (FEA) for elasticity / heat conduction,

• density initialization and masking of fixed / boundary elements,

• density and sensitivity filtering, and

• Heaviside-type projection and sensitivity chain rules,

while leaving the actual density update rule to subclasses via rho_update().

Typical usage involves:

• Assembling the global stiffness / conductivity matrix (via an internal FEM_* instance),

• Solving the state problem under given loads and boundary conditions,

• Computing compliance (or thermal) objectives and sensitivities,

• Applying density / sensitivity filters (e.g. spatial, Helmholtz),

• Projecting intermediate densities with a smooth Heaviside function, and

• Delegating the update of design densities to a subclass-specific implementation (OC, MOC variants, etc.).

Responsibilities

• Manage the overall optimization loop (iterations, restart, output).

• Initialize and update continuation schedules (penalization, volume fraction, projection sharpness, move limits, etc.).

• Apply filtering and projection consistently to densities and gradients.

• Evaluate objective values and record histories / visualization outputs.

• Call rho_update() to perform the actual density update (including, e.g., Lagrange multiplier bisection for OC).

Notes

This class serves as the backbone of density-based topology optimization in Scikit-Topt. All algorithm-specific update logic (e.g., OC / MOC density updates, dual variable handling) must be implemented in subclasses by overriding rho_update().

cfg

Configuration object holding numerical and algorithmic settings, such as interpolation model, filters, continuation schedules, and solver options.

Type:

DensityMethodConfig

tsk

Problem description including mesh, basis, boundary conditions, and load vectors.

Type:

sktopt.mesh.FEMDomain

fem

Finite element solver used for state and sensitivity evaluation, constructed from tsk and cfg.

Type:

sktopt.fea.FEM_SimpLinearElasticity or sktopt.fea.FEM_SimpLinearHeatConduction

schedulers

Collection of SchedulerConfig instances controlling continuation of parameters such as p, vol_frac, beta, move_limit, percentile, and filter_radius.

Type:

sktopt.tools.Schedulers

initialize_params()

Initialize parameters and preallocate work arrays for density-based topology optimization using the SIMP method.

This prepares all arrays required for the iterative loop. By preallocating and reusing arrays in-place, repeated allocations are avoided and performance improves.

Returns:

• iter_begin (int) – Starting iteration index.

• iter_end (int) – Final iteration index.

• rho (ndarray, shape (n_elements,)) – Current element-wise design density.

• rho_prev (ndarray, shape (n_elements,)) – Copy of densities from the previous iteration.

• rho_filtered (ndarray, shape (n_elements,)) – Densities after applying the density filter (e.g., Helmholtz filter).

• rho_projected (ndarray, shape (n_elements,)) – Densities after applying the projection (e.g., Heaviside projection).

• dH_drho (ndarray, shape (n_elements,)) – Derivative of the projection with respect to the filtered densities.

• grad_filtered (ndarray, shape (n_elements,)) – Aggregated gradient after combining sensitivity with projection derivative.

• dC_drho_projected (ndarray, shape (n_elements,)) – Compliance sensitivities with respect to the projected densities.

• energy_mean (ndarray, shape (n_elements,)) – Average element strain energy over all load cases.

• dC_drho_full (ndarray, shape (n_elements,)) – Compliance sensitivities mapped back to the full element set.

• dC_drho_design_eles (ndarray, shape (n_design_elements,)) – Compliance sensitivities restricted to design elements only.

• scaling_rate (ndarray, shape (n_design_elements,)) – Scaling factors used in the update scheme (e.g., OC/MOC).

• rho_design_eles (ndarray, shape (n_design_elements,)) – Subset of densities corresponding to design elements only.

• rho_clip_lower (ndarray, shape (n_design_elements,)) – Lower clipping bounds for density updates.

• rho_clip_upper (ndarray, shape (n_design_elements,)) – Upper clipping bounds for density updates.

• u_dofs (ndarray, shape (ndof, n_load_cases)) – Displacement field solutions for each degree of freedom and each load case.

• filter_radius (float) – Initial density-filter radius, determined from a fixed value or a schedule.

Notes

• This setup targets density-based (SIMP) topology optimization.

• Arrays are reused in-place across iterations to minimize allocation overhead.

optimize()

Run optimization until cfg.max_iters (default behavior).

optimize_steps(num_steps: int)

Run only num_steps iterations of the optimizer and record progress.

This is useful when an external driver wants to probe or checkpoint intermediate states without waiting for the full max_iters loop.

class sktopt.core.DensityMethodConfig(dst_path: str = './result/pytests', interpolation: ~typing.Literal['SIMP', 'RAMP'] = 'SIMP', record_times: int = 20, max_iters: int = 200, beta_eta: float = 0.5, eta: float = 0.6, p: ~sktopt.tools.SchedulerConfig = <factory>, vol_frac: ~sktopt.tools.SchedulerConfig = <factory>, beta: ~sktopt.tools.SchedulerConfig = <factory>, neumann_scale: ~sktopt.tools.SchedulerConfig = <factory>, filter_type: ~typing.Literal['spacial', 'helmholtz'] = 'helmholtz', filter_radius: ~sktopt.tools.SchedulerConfig = <factory>, E_min_coeff: float = 0.001, rho_min: float = 0.01, rho_max: float = 1.0, restart: bool = False, restart_from: int = -1, export_img: bool = False, export_img_opaque: bool = False, design_dirichlet: bool = False, sensitivity_filter: bool = False, solver_option: ~typing.Literal['spsolve', 'cg_pyamg'] = 'spsolve', scaling: bool = False, check_convergence: bool = False, tol_rho_change: float = 0.2, tol_kkt_residual: float = 0.005)

Bases: object

Configuration for density-based topology optimization (SIMP).

This configuration collects the main numerical settings for density-based topology optimization, including interpolation models, continuation schedules (via SchedulerConfig), filtering/projection parameters, and solver options. It is intended for sensitivity-based methods such as SIMP and RAMP and is shared by both OC and MOC-style update rules.

Notes on continuation

Parameters that change during optimization are represented by SchedulerConfig:

• p: Penalization power (e.g., 1.0 → 3.0).

• vol_frac: Target volume fraction (often constant, e.g. 0.8).

• beta: Heaviside / projection sharpness (e.g., 1.0 → 2.0 or higher).

• filter_radius: Filter radius (typically kept constant).

Each scheduler defines how its value evolves during the optimization (e.g., SchedulerConfig.step(), SchedulerConfig.step_accelerating(), SchedulerConfig.constant()) and may include shape parameters such as curvature for accelerating schedules.

dst_path

Output directory for results and logs.

Type:

str

interpolation

Material interpolation model used for penalization.

Type:

{“SIMP”, “RAMP”}

record_times

Number of snapshots recorded during optimization.

Type:

int

max_iters

Maximum number of optimization iterations.

Type:

int

beta_eta

Additional damping parameter used in OC/MOC-style updates (for example, to smooth dual variables or stabilize multiplicative/log-space updates).

Type:

float

eta

Damping or learning-rate–like exponent used in the OC/MOC multiplicative update rule.

Type:

float

p

Continuation schedule for the penalization power. Default: step(init_value=1.0, target_value=3.0, num_steps=3).

Type:

SchedulerConfig

vol_frac

Schedule for the target volume fraction constraint. Default: constant(target_value=0.8).

Type:

SchedulerConfig

beta

Continuation schedule for Heaviside / projection sharpness (β). Default: step_accelerating(init_value=1.0, target_value=2.0, num_steps=3, curvature=2.0).

Type:

SchedulerConfig

neumann_scale

Continuation scale applied multiplicatively to the Neumann force vector. Only ConstantOne or StepToOne / Step*ToOne are allowed. Default: constant 1.0 (no scaling).

Type:

SchedulerConfig

filter_type

Type of filter applied to the density field (spatial or Helmholtz).

Type:

{“spacial”, “helmholtz”}

filter_radius

Schedule for the filter radius. Default: constant(target_value=0.01).

Type:

SchedulerConfig

E_min_coeff

Proportional constant that defines the minimum Young’s modulus as E_min_coeff * E0 (typically \(10^{-3}\)).

Type:

float

rho_min

Minimum density clamp to avoid singular stiffness matrices.

Type:

float

rho_max

Maximum density clamp (typically 1.0).

Type:

float

restart

If True, resume optimization from a saved state in dst_path.

Type:

bool

restart_from

Iteration index to resume from. Use -1 to auto-detect the latest checkpoint.

Type:

int

export_img

If True, export images of the density field during optimization.

Type:

bool

export_img_opaque

If True, use opaque rendering for exported images.

Type:

bool

design_dirichlet

If True, Dirichlet boundary elements are included in the design domain.

Type:

bool

sensitivity_filter

If True, apply filtering directly to sensitivity fields in addition to density filtering.

Type:

bool

solver_option

Linear solver for the state analysis. "cg_pyamg" enables multigrid-accelerated CG via PyAMG.

Type:

{“spsolve”, “cg_pyamg”}

scaling

If True, apply length/force scaling to normalize geometry and loads for improved numerical conditioning.

Type:

bool

check_convergence

If True, enable automatic convergence checking during optimization. Convergence is determined by two criteria:

1. Maximum density change (tol_rho_change)

2. KKT-like residual or Lagrangian gradient norm (tol_kkt_residual)

Optimization terminates early only when both are satisfied.

Type:

bool

tol_rho_change

Tolerance for the maximum per-iteration density change. Convergence is considered achieved when:

max(|rho_new - rho_old|) < tol_rho_change

Typical values range from 1e-2 to 1e-4 depending on mesh resolution.

Type:

float

tol_kkt_residual

Tolerance for the KKT-style stationarity residual.

For OC:

The residual is computed as the infinity norm of the Lagrangian gradient dC/drho + λ * dV/drho over interior design elements.

For MOC / LogMOC:

A Lagrangian-like gradient norm based on the current volume-penalty multiplier is used instead, serving as an approximate stationarity measure.

Convergence is considered achieved when:

kkt_residual < tol_kkt_residual

Type:

float

class sktopt.core.DensityMethod_OC_Config(dst_path: str = './result/pytests', interpolation: ~typing.Literal['SIMP', 'RAMP'] = 'SIMP', record_times: int = 20, max_iters: int = 200, beta_eta: float = 0.5, eta: float = 0.6, p: ~sktopt.tools.SchedulerConfig = <factory>, vol_frac: ~sktopt.tools.SchedulerConfig = <factory>, beta: ~sktopt.tools.SchedulerConfig = <factory>, neumann_scale: ~sktopt.tools.SchedulerConfig = <factory>, filter_type: ~typing.Literal['spacial', 'helmholtz'] = 'helmholtz', filter_radius: ~sktopt.tools.SchedulerConfig = <factory>, E_min_coeff: float = 0.001, rho_min: float = 0.01, rho_max: float = 1.0, restart: bool = False, restart_from: int = -1, export_img: bool = False, export_img_opaque: bool = False, design_dirichlet: bool = False, sensitivity_filter: bool = False, solver_option: ~typing.Literal['spsolve', 'cg_pyamg'] = 'spsolve', scaling: bool = False, check_convergence: bool = False, tol_rho_change: float = 0.2, tol_kkt_residual: float = 0.005, lambda_lower: float = 1e-07, lambda_upper: float = 10000000.0, percentile: ~sktopt.tools.SchedulerConfig = <factory>, move_limit: ~sktopt.tools.SchedulerConfig = <factory>)

Bases: DensityMethodConfig

Configuration class for OC-style density updates.

This subclass of DensityMethodConfig adds parameters that are specific to Optimality Criteria (OC) and related multiplicative update schemes. It controls the move-limit strategy, optional percentile-based sensitivity scaling, and the bracket for the Lagrange multiplier used in the volume-constraint bisection.

lambda_lower

Lower bound for the Lagrange multiplier used in the bisection enforcing the volume constraint. Must be strictly positive.

Type:

float

lambda_upper

Upper bound for the Lagrange multiplier used in the bisection enforcing the volume constraint.

Type:

float

percentile

Optional schedule for percentile-based scaling of sensitivity fields. When set to SchedulerConfig.none(), percentile scaling is disabled and raw sensitivities are used.

Type:

SchedulerConfig

move_limit

Schedule for the OC move limit applied to density updates. By default, a sawtooth-decay schedule is used, e.g. sawtooth_decay("move_limit", 0.3, 0.1, 3), which gradually tightens the allowed per-iteration change in density.

Type:

SchedulerConfig

class sktopt.core.DensityState(rho: numpy.ndarray, rho_prev: numpy.ndarray, rho_filtered: numpy.ndarray, rho_projected: numpy.ndarray, dH_drho: numpy.ndarray, grad_filtered: numpy.ndarray, dC_drho_projected: numpy.ndarray, energy_mean: numpy.ndarray, dC_drho_full: numpy.ndarray, dC_drho_design_eles: numpy.ndarray, scaling_rate: numpy.ndarray, rho_design_eles: numpy.ndarray, rho_clip_lower: numpy.ndarray, rho_clip_upper: numpy.ndarray, u_dofs: numpy.ndarray, filter_radius: float, elements_volume_design: numpy.ndarray, elements_volume_design_sum: float, iter_begin: int, iter_end: int, last_iter: int | None = None, compliance: float | None = None, u_max: float | numpy.ndarray | None = None, rho_change_max: float | None = None, kkt_residual: float | None = None, vol_error: float | None = None)

Bases: object

class sktopt.core.LogMOC_Config(dst_path: str = './result/pytests', interpolation: ~typing.Literal['SIMP'] = 'SIMP', record_times: int = 20, max_iters: int = 200, beta_eta: float = 0.5, eta: float = 0.6, p: ~sktopt.tools.SchedulerConfig = <factory>, vol_frac: ~sktopt.tools.SchedulerConfig = <factory>, beta: ~sktopt.tools.SchedulerConfig = <factory>, neumann_scale: ~sktopt.tools.SchedulerConfig = <factory>, filter_type: ~typing.Literal['spacial', 'helmholtz'] = 'helmholtz', filter_radius: ~sktopt.tools.SchedulerConfig = <factory>, E_min_coeff: float = 0.001, rho_min: float = 0.01, rho_max: float = 1.0, restart: bool = False, restart_from: int = -1, export_img: bool = False, export_img_opaque: bool = False, design_dirichlet: bool = False, sensitivity_filter: bool = False, solver_option: ~typing.Literal['spsolve', 'cg_pyamg'] = 'spsolve', scaling: bool = False, check_convergence: bool = False, tol_rho_change: float = 0.2, tol_kkt_residual: float = 0.005, lambda_lower: float = 1e-07, lambda_upper: float = 10000000.0, percentile: ~sktopt.tools.SchedulerConfig = <factory>, move_limit: ~sktopt.tools.SchedulerConfig = <factory>, mu_p: float = 0.1, lambda_v: float = 100.0, lambda_decay: float = 0.7, scaling_rate_min: float = -0.5, scaling_rate_max: float = 0.5)

Bases: DensityMethod_OC_Config

Configuration for a log-space, OC-like multiplicative density update.

This configuration supports a variant of the Optimality Criteria (OC) method implemented in \(\log \rho\)-space. The update remains multiplicative in \(\rho\), but it is applied by adding a scaled log-factor to \(\log \rho\). The dual variable associated with the volume constraint is updated by an exponential moving average (EMA) driven by the volume error.

Conceptually, the update has the form

\[\rho^{\text{new}} = \rho^{\text{old}} \left(-\frac{\partial C / \partial \rho}{\lambda_v} \right)^{\eta},\]

where the ratio \(-\partial C/\partial \rho / \lambda_v\) is clamped to avoid numerical instabilities and move limits are enforced in log-space.

interpolation

Interpolation scheme used for material penalization. Currently fixed to "SIMP".

Type:

Literal[“SIMP”]

mu_p

Gain used to convert the volume error into a dual update. Larger values make the dual variable lambda_v react more strongly to constraint violations.

Type:

float

lambda_v

Initial value for the volume-related dual variable. This is used in the OC-like scaling term \(-\partial C/\partial \rho / \lambda_v\).

Type:

float

lambda_decay

Decay factor in the EMA update of lambda_v. Values closer to 1.0 yield slower, smoother updates; smaller values react more quickly to changes in the volume error.

Type:

float

class sktopt.core.LogMOC_Optimizer(cfg: LogMOC_Config, tsk: FEMDomain)

Bases: DensityMethod

Topology optimization solver using a log-space, OC-like multiplicative update.

This optimizer implements a density-based compliance minimization algorithm in which the classic OC multiplicative update is expressed in log(ρ)-space. The design densities are updated according to a smoothed ratio of sensitivities and a dual variable associated with the volume constraint.

In terms of ρ, the update has the conceptual form

\[\rho^{\text{new}} = \rho^{\text{old}} \left(-\frac{\partial C / \partial \rho}{\lambda_v} \right)^{\eta},\]

where \(\lambda_v\) is a dual variable updated by an exponential moving average (EMA) driven by the volume constraint violation. The implementation performs this update in log-space, enforces move limits on \(\log \rho\), and finally clips the densities to the admissible interval \([\rho_{\min}, \rho_{\max}]\).

Compared to the classic OC method, this log-space variant can offer improved numerical robustness while preserving the intuitive multiplicative structure.

cfg

Configuration object specifying the OC-like settings, including mu_p, lambda_v, EMA decay factors, and bounds on the densities and dual variables.

Type:

LogMOC_Config

tsk

Finite element model (mesh, basis, boundary conditions and loads).

Type:

sktopt.mesh.FEMDomain

recorder

Logger that stores histories of selected quantities such as -dC, lambda_v, volume error and KKT residuals.

Type:

HistoryLogger

class sktopt.core.OC_Config(dst_path: str = './result/pytests', interpolation: ~typing.Literal['SIMP'] = 'SIMP', record_times: int = 20, max_iters: int = 200, beta_eta: float = 0.5, eta: ~sktopt.tools.SchedulerConfig = <factory>, p: ~sktopt.tools.SchedulerConfig = <factory>, vol_frac: ~sktopt.tools.SchedulerConfig = <factory>, beta: ~sktopt.tools.SchedulerConfig = <factory>, neumann_scale: ~sktopt.tools.SchedulerConfig = <factory>, filter_type: ~typing.Literal['spacial', 'helmholtz'] = 'helmholtz', filter_radius: ~sktopt.tools.SchedulerConfig = <factory>, E_min_coeff: float = 0.001, rho_min: float = 0.01, rho_max: float = 1.0, restart: bool = False, restart_from: int = -1, export_img: bool = False, export_img_opaque: bool = False, design_dirichlet: bool = False, sensitivity_filter: bool = False, solver_option: ~typing.Literal['spsolve', 'cg_pyamg'] = 'spsolve', scaling: bool = False, check_convergence: bool = False, tol_rho_change: float = 0.2, tol_kkt_residual: float = 0.005, lambda_lower: float = 1e-07, lambda_upper: float = 10000000.0, percentile: ~sktopt.tools.SchedulerConfig = <factory>, move_limit: ~sktopt.tools.SchedulerConfig = <factory>, scaling_rate_min: float = 0.7, scaling_rate_max: float = 1.3)

Bases: DensityMethod_OC_Config

Optimality-Criteria (OC) configuration for density-based optimization.

This class specializes common_density.DensityMethod_OC_Config for OC-style updates. The material interpolation is fixed to SIMP, and the projection threshold eta is controlled via a sktopt.tools.SchedulerConfig.

interpolation

Material interpolation model. Fixed to “SIMP” for OC.

Type:

{“SIMP”}

eta

Threshold (η) used in Heaviside projection / update logic. The default keeps η constant at 0.5 using a Constant scheduler. Pass a custom scheduler to run continuation (e.g., Step from 0.7 → 0.3).

Type:

sktopt.tools.SchedulerConfig

Notes

• Other continuation parameters (e.g., p, beta, vol_frac, filter_radius) are inherited from the parent config if present there.

• To override eta from outside, provide a pre-built SchedulerConfig at construction time.

class sktopt.core.OC_Optimizer(cfg: OC_Config, tsk: FEMDomain)

Bases: DensityMethod

Topology optimization solver using the classic Optimality Criteria (OC) method. This class implements the standard OC algorithm for compliance minimization problems. It uses a multiplicative density update formula derived from Karush-Kuhn-Tucker (KKT) optimality conditions under volume constraints.

The update rule typically takes the form:

ρ_new = clamp(ρ * sqrt(-dC / λ), ρ_min, ρ_max)

where:

• dC is the sensitivity of the compliance objective,

• λ is a Lagrange multiplier for the volume constraint.

This method is widely used in structural optimization due to its simplicity, interpretability, and solid theoretical foundation.

Advantages

• Simple and easy to implement

• Intuitive update rule based on physical insight

• Well-established and widely validated in literature

config

Configuration object specifying the interpolation method, volume fraction, continuation settings, filter radius, and other numerical parameters.

Type:

DensityMethodConfig

mesh, basis, etc.

FEM components required for simulation, including boundary conditions and loads.

Type:

inherited from common_density.DensityMethod