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', 'petsc', 'petsc_spdirect'] = 'spsolve', petsc_options: ~sktopt.fea.solver_petsc.PETScOptions = <factory>, 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

Legacy convenience selector for the state-analysis solver. "cg_pyamg" enables multigrid-accelerated CG via PyAMG. "petsc" uses PETSc CG + GAMG (requires petsc4py). "petsc_spdirect" uses PETSc sparse direct solve (ksp_type=preonly, pc_type=lu).

Type:: {“spsolve”, “cg_pyamg”, “petsc”, “petsc_spdirect”}

petsc_options

Legacy convenience PETSc options used to build solver_config. Additional PETSc-specific options are forwarded only when solver_option is "petsc" or "petsc_spdirect". Supported keys currently include ksp_type, pc_type, pc_factor_mat_solver_type, options_prefix, and use_set_from_options.

Type:: PETScOptions

solver_config

Normalized solver configuration derived from solver_option and petsc_options during initialization.

Type:: LinearSolverConfig

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:

Maximum density change (tol_rho_change)
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', 'petsc', 'petsc_spdirect'] = 'spsolve', petsc_options: ~sktopt.fea.solver_petsc.PETScOptions = <factory>, 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

sktopt.core.LogMOCScaled_Config: alias of LogMOC_Config

sktopt.core.LogMOCScaled_Optimizer: alias of LogMOC_Optimizer

class sktopt.core.LogMOC_Config(dst_path: str = './result/pytests', interpolation: 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: 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: Literal['spsolve', 'cg_pyamg', 'petsc', 'petsc_spdirect'] = 'spsolve', petsc_options: sktopt.fea.solver_petsc.PETScOptions = <factory>, scaling: bool = False, check_convergence: bool = False, tol_rho_change: float = 0.2, tol_kkt_residual: float = 0.005, lambda_lower: float = -10000000.0, lambda_upper: float = 10000000.0, percentile: sktopt.tools.SchedulerConfig = <factory>, move_limit: sktopt.tools.SchedulerConfig = <factory>, mu_p: float = 5.0, augmented_lagrangian_mu: float = 0.0, lambda_v: float = 0.1, lambda_decay: float = 0.9, lagrangian_clip: float = 1.0, lagrangian_percentile: float = 95.0, lagrangian_scale_floor: float = 1e-08, normalize_volume_chain: bool = False, volume_chain_percentile: float = 95.0, volume_chain_scale_floor: float = 1e-08, volume_chain_gain: float = 1.0, volume_chain_gain_under: float | None = None, dual_update: Literal['ema', 'augmented'] = 'ema', filter_lagrangian: bool = False, center_lagrangian: bool = False, center_objective: bool = False): Bases: DensityMethod_OC_Config

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

sktopt.core.OCScaled_Config: alias of OC_Config

sktopt.core.OCScaled_Optimizer: alias of OC_Optimizer

class sktopt.core.OC_Config(dst_path: str = './result/pytests', interpolation: 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: 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: Literal['spsolve', 'cg_pyamg', 'petsc', 'petsc_spdirect'] = 'spsolve', petsc_options: sktopt.fea.solver_petsc.PETScOptions = <factory>, 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, sensitivity_scale_floor: float = 1e-12): Bases: DensityMethod_OC_Config

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

Bases: DensityMethod

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.