API Reference

This section provides the API documentation for the core modules in the ggp framework.

Problem Specification

Serialisable dataclass specifications for GGP optimisation problems.

Every optimisation run is fully described by a single ProblemSpec instance. This object can be constructed programmatically, loaded from a YAML file, or built from CLI arguments.

class ggp.problem.spec.BoundaryCondition(region: str, type: Literal['fixed', 'symmetry', 'prescribed'] = 'fixed', components: List[int] | None = None, value: List[float] | None = None)[source]

Bases: object

A Dirichlet boundary condition applied to a named region.

Parameters

regionstr: Boundary identifier (e.g. "left", "top", "fixed_face"). Mapped to concrete selection logic by the discretiser.
type{“fixed”, “symmetry”, “prescribed”}: Kind of constraint.
componentslist of int, optional: DOF components to constrain (0 = x, 1 = y, 2 = z). None means all components.
valuelist of float, optional: Prescribed displacement for type="prescribed".

components: List[int] | None = None

region: str

type: Literal['fixed', 'symmetry', 'prescribed'] = 'fixed'

value: List[float] | None = None

class ggp.problem.spec.ConstraintSpec(name: str, type: ~typing.Literal['ineq', 'eq'] = 'ineq', bound: float = 0.0, params: ~typing.Dict[str, ~typing.Any] = <factory>)[source]

Bases: object

An optimisation constraint.

Parameters

namestr: Constraint identifier ("volume", "overhang", "stress", …).
type{“ineq”, “eq”}: Inequality or equality.
boundfloat: Right-hand-side value (constraint ≤ bound for ineq).
paramsdict: Extra parameters (e.g. {"alpha_deg": 45} for overhang).

bound: float = 0.0

name: str

params: Dict[str, Any]

type: Literal['ineq', 'eq'] = 'ineq'

class ggp.problem.spec.FormulationSpec(mode: Literal['Free', 'ALM', '3D_Free', '3D_ALM', 'Truss', '2D_Truss', 'Truss3D', '3D_Truss'] = 'Free', num_components: int = 18, num_layers: int | None = None, comp_per_layer: int | None = None, layer_height: float | None = None, r_gp: float = 0.5, Ngp: int = 2, method: str = 'GP', init: str = 'default', grid_nx: int | None = None, grid_ny: int | None = None, min_thickness: float | None = None, symmetry: str | None = None, truss_radius: float | None = None, grid_nz: int | None = None, fix_mc: bool = False)[source]

Bases: object

GGP projection formulation specification.

Parameters

mode{“Free”, “ALM”, “3D_Free”, “3D_ALM”}: Projection mode.
num_componentsint: Total number of geometric primitives.
num_layersint, optional: Number of manufacturing layers (ALM modes only).
comp_per_layerint, optional: Components per layer (ALM modes only).
layer_heightfloat, optional: Height of each layer (ALM modes only).
r_gpfloat: Sampling-window radius for Gauss-point integration.
Ngpint: Number of Gauss points per dimension.
method{“GP”, “WGP”, “rect”}: Geometry projection method / primitive variant. "GP" is the round-ended bar (capsule); "rect" is a true rectangle primitive (quintic-smoothed, width L x height h).
init{“default”, “grid”}: Initial-design layout. "default" is the Matlab-style crossed-bar grid; "grid" tiles the domain with grid_nx x grid_ny rectangles (squares when the tile aspect matches the domain) at density = volfrac.
grid_nx, grid_nyint, optional: Number of tiles along x and y for the "grid" initialisation. Their product must equal num_components and both must be < the mesh resolution.

Ngp: int = 2

comp_per_layer: int | None = None

fix_mc: bool = False

grid_nx: int | None = None

grid_ny: int | None = None

grid_nz: int | None = None

init: str = 'default'

layer_height: float | None = None

method: str = 'GP'

min_thickness: float | None = None

mode: Literal['Free', 'ALM', '3D_Free', '3D_ALM', 'Truss', '2D_Truss', 'Truss3D', '3D_Truss'] = 'Free'

num_components: int = 18

num_layers: int | None = None

r_gp: float = 0.5

symmetry: str | None = None

truss_radius: float | None = None

class ggp.problem.spec.GeometrySpec(type: str, role: ~typing.Literal['design', 'non_design'] = 'design', path: ~pathlib.Path | None = None, params: ~typing.Dict[str, ~typing.Any] = <factory>)[source]

Bases: object

A geometry region (design-space or non-design-space).

Supported type values (extensible via the reader registry):

type	Description
`box`	Parametric axis-aligned box (2-D or 3-D)
`fenics_rectangle`	FEniCS built-in quadrilateral rectangle
`fenics_box`	FEniCS built-in hexahedral box
`step`	STEP CAD file (requires `gmsh` or `cadquery`)
`stl`	STL surface mesh (requires `meshio`)
`brep`	OpenCASCADE BREP file
`points`	Raw point cloud (`numpy` only)

Parameters

typestr: Geometry source type.
role{“design”, “non_design”}: Whether this region participates in the design optimisation.
pathPath, optional: File path for file-based formats (STEP, STL, …).
paramsdict: Format-specific parameters (e.g. {"Lx": 60, "Ly": 30, "nx": 60, "ny": 30}).

params: Dict[str, Any]

path: Path | None = None

role: Literal['design', 'non_design'] = 'design'

type: str

class ggp.problem.spec.Load(region: str, type: ~typing.Literal['point', 'patch', 'pressure', 'traction', 'body'] = 'point', value: ~typing.List[float] = <factory>, width: float = 0.0)[source]

Bases: object

A mechanical load applied to a named region.

Parameters

regionstr: Region identifier (e.g. "mid_right", "top_face").
type{“point”, “patch”, “pressure”, “traction”, “body”}: Load application mode. "patch" spreads the total force uniformly over the boundary nodes within width of the region’s reference point — the standard way to avoid the point-load stress singularity in stress-constrained problems.
valuelist of float: Force/pressure vector (total force for "point"/"patch").
widthfloat: Extent of the "patch" load along the loaded face (ignored otherwise).

region: str

type: Literal['point', 'patch', 'pressure', 'traction', 'body'] = 'point'

value: List[float]

width: float = 0.0

class ggp.problem.spec.MaterialSpec(E: float = 1.0, nu: float = 0.3)[source]

Bases: object

Isotropic material properties.

Parameters

Efloat: Young’s modulus.
nufloat: Poisson’s ratio.

E: float = 1.0

nu: float = 0.3

class ggp.problem.spec.ObjectiveSpec(name: str = 'compliance', sense: ~typing.Literal['minimize', 'maximize'] = 'minimize', params: ~typing.Dict[str, ~typing.Any] = <factory>)[source]

Bases: object

Objective function to minimise or maximise.

Parameters

namestr: Response identifier — one of KNOWN_RESPONSES ("compliance", "volume", "displacement", "stress").
sense{“minimize”, “maximize”}: Optimisation direction.
paramsdict: Extra parameters for responses that need them when used as the objective without also being a constraint (e.g. P/aggregation for "stress", axis/point for "displacement") — mirrors ConstraintSpec.params.

name: str = 'compliance'

params: Dict[str, Any]

sense: Literal['minimize', 'maximize'] = 'minimize'

class ggp.problem.spec.ProblemSpec(geometries: ~typing.List[~ggp.problem.spec.GeometrySpec], boundary_conditions: ~typing.List[~ggp.problem.spec.BoundaryCondition], loads: ~typing.List[~ggp.problem.spec.Load], formulation: ~ggp.problem.spec.FormulationSpec, objective: ~ggp.problem.spec.ObjectiveSpec = <factory>, constraints: ~typing.List[~ggp.problem.spec.ConstraintSpec] = <factory>, solver: ~ggp.problem.spec.SolverSpec = <factory>, volfrac: float = 0.4, material: ~ggp.problem.spec.MaterialSpec = <factory>)[source]

Bases: object

Complete, self-contained optimisation problem definition.

This is the single source of truth for any GGP run. The CLI, Python API, and YAML loader all produce a ProblemSpec.

Parameters

geometrieslist of GeometrySpec: Domain geometry regions (design + non-design).
boundary_conditionslist of BoundaryCondition: Dirichlet boundary conditions.
loadslist of Load: Applied loads.
formulationFormulationSpec: GGP formulation / projection configuration.
objectiveObjectiveSpec: Objective function.
constraintslist of ConstraintSpec: Optimisation constraints.
solverSolverSpec: Algorithm and solver configuration.
volfracfloat: Target volume fraction.
materialMaterialSpec: Material properties.

boundary_conditions: List[BoundaryCondition]

constraints: List[ConstraintSpec]

formulation: FormulationSpec

geometries: List[GeometrySpec]

loads: List[Load]

material: MaterialSpec

objective: ObjectiveSpec

solver: SolverSpec

volfrac: float = 0.4

class ggp.problem.spec.SolverSpec(algorithm: str = 'MMA', max_iter: int = 50, iterative: bool = False, fem_solver: ~typing.Literal['direct', 'iterative', 'amjax'] = 'direct', use_line_search: bool = False, options: ~typing.Dict[str, ~typing.Any] = <factory>)[source]

Bases: object

Optimisation solver / algorithm configuration.

Parameters

algorithmstr: Optimisation algorithm name ("MMA", "SLP", "CONLIN").
max_iterint: Maximum number of outer iterations.
iterativebool: Deprecated shorthand for fem_solver="iterative". Kept for backward compatibility with existing YAML files and CLI.
fem_solver{“direct”, “iterative”, “amjax”}: FEM linear-system backend. "direct" uses scipy.sparse.linalg.spsolve() (LU); "iterative" uses PETSc CG + GAMG; "amjax" uses PyAMG + AMJax (JAX-accelerated algebraic multigrid, JIT-compiled and GPU-compatible).
use_line_searchbool: Enable line search (SLP / CONLIN only).
optionsdict: Additional algorithm-specific key-value options.

algorithm: str = 'MMA'

fem_solver: Literal['direct', 'iterative', 'amjax'] = 'direct'

iterative: bool = False

max_iter: int = 50

options: Dict[str, Any]

use_line_search: bool = False

Geometry I/O & Parsers

Self-registering geometry reader registry.

Usage

Register a reader:

@register_reader("step")
class STEPReader(GeometryReader):
    ...

Retrieve a reader:

reader = get_reader("step")
domain = reader.read(geometry_spec)

List available readers:

names = list_readers()  # ["fenics_rectangle", "fenics_box", "step", ...]

ggp.geometry.io.registry.get_reader(geometry_type: str) → GeometryReader[source]

Instantiate the registered reader for geometry_type.

Raises

ValueError: If geometry_type has not been registered.

ggp.geometry.io.registry.list_readers() → List[str][source]: Return the names of all registered geometry readers.

ggp.geometry.io.registry.register_reader(geometry_type: str)[source]

Class decorator that registers a GeometryReader subclass.

Parameters

geometry_typestr: The GeometrySpec.type value this reader handles.

FEniCS built-in geometry reader.

Handles fenics_rectangle (2-D quadrilateral mesh) and fenics_box (3-D hexahedral mesh) geometry types by calling FEniCS mesh generators.

class ggp.geometry.io.fenics_reader.FenicsBoxReader[source]

Bases: GeometryReader

Creates a structured hexahedral box via dolfin.BoxMesh.

read(spec: Any) → DomainRepresentation[source]

Convert spec (a GeometrySpec) into a domain representation.

Parameters

specGeometrySpec: Geometry specification from the problem definition layer.

Returns

DomainRepresentation

static supports(geometry_type: str) → bool[source]: Return True if this reader handles geometry_type.

class ggp.geometry.io.fenics_reader.FenicsRectangleReader[source]

Bases: GeometryReader

Creates a structured quadrilateral rectangle via dolfin.RectangleMesh.

read(spec: Any) → DomainRepresentation[source]

Convert spec (a GeometrySpec) into a domain representation.

Parameters

specGeometrySpec: Geometry specification from the problem definition layer.

Returns

DomainRepresentation

static supports(geometry_type: str) → bool[source]: Return True if this reader handles geometry_type.

Discretisation

Discretisation layer — FEM mesh discretiser.

class ggp.discretisation.fem.AnalysisDomain(dim: int, eval_coords: numpy.ndarray, metadata: dict | None = None)[source]

Bases: object

Analysis-ready domain containing mesh and BCs.

class ggp.discretisation.fem.FEMDiscretiser[source]

Bases: object

Converts a DomainRepresentation into a FEniCS-ready AnalysisDomain.

discretise(domain: DomainRepresentation, spec: ProblemSpec) → AnalysisDomain[source]

Projection Mappers

Self-registering mapper registry.

Usage

Register a mapper:

@register_mapper("Free")
class Free2DMapper(ProjectionMapper):
    ...

Retrieve a mapper:

mapper = get_mapper("Free", num_components=18, ...)

List available mappers:

names = list_mappers()  # ["Free", "ALM", "3D_Free", ...]

ggp.projection.registry.get_mapper(mode: str, **kwargs) → ProjectionMapper[source]

Instantiate the registered mapper for mode.

Parameters

modestr: Formulation mode (must match a previously registered mapper).
**kwargs: Forwarded to the mapper constructor.

Raises

ValueError: If mode has not been registered.

ggp.projection.registry.list_mappers() → List[str][source]: Return the names of all registered mappers.

ggp.projection.registry.register_mapper(mode: str)[source]

Class decorator that registers a ProjectionMapper subclass.

Parameters

modestr: Formulation mode identifier (e.g. "Free", "ALM", "3D_Free").

Free 2-D geometry projection mapper.

Maps 6 variables per component: [Xc, Yc, L, h, Theta, Mc].

class ggp.projection.free_2d.Free2DMapper(num_components: int, r_gp: float = 0.5, method: str = 'GP', Ngp: int = 2, min_thickness: float = 1.0, fix_mc: bool = False, **kwargs)[source]

Bases: ProjectionMapper

Free 2-D projection mapper.

Each component uses 6 parameters:: [Xc, Yc, length, thickness, angle, density]

default_bounds(domain_extents: Tuple[float, ...], num_components: int) → Tuple[numpy.ndarray, numpy.ndarray][source]

Return (lb, ub) arrays for the full design-variable vector.

Parameters

domain_extentstuple of float: (Lx, Ly) for 2-D, (Lx, Ly, Lz) for 3-D.
num_componentsint: Number of geometric primitives.

Returns

lb, ubndarray, ndarray: Lower and upper bounds, each of shape (n_vars,).

forward(x_vars: numpy.ndarray, eval_coords: numpy.ndarray, power_E: float = 1.0, power_V: float = 1.0) → Tuple[numpy.ndarray, numpy.ndarray][source]

Compute density fields.

Parameters

x_varsndarray, shape (n_vars,): Flat design-variable vector (already unscaled).
eval_coordsndarray, shape (n_eval, dim): Element centroid (or evaluation-point) coordinates.
power_Efloat: Penalisation exponent for stiffness density.
power_Vfloat: Penalisation exponent for volume density.

Returns

rho_Endarray, shape (n_elements,): Stiffness density per element.
rho_Vndarray, shape (n_elements,): Volume density per element.

jacobian(x_vars: numpy.ndarray, eval_coords: numpy.ndarray, power_E: float = 1.0, power_V: float = 1.0) → Dict[str, numpy.ndarray][source]

Compute analytic Jacobians.

Returns

dict: {"rho_E": drho_E_dx, "rho_V": drho_V_dx} where each value has shape (n_elements, n_vars).

num_vars_per_component() → int[source]: Number of design variables per geometric primitive.

Additive Layer Manufacturing (ALM) 2-D geometry projection mapper.

Correct variable layout (interleaved, matching GGP_ALM_Matlab reference):

[Xc_0_0, L_0_0, Xc_0_1, L_0_1, …, Xc_{nY-1}_{np-1}, L_{nY-1}_{np-1},: h_0, …, h_{np-1}, Mc_0, …, Mc_{np-1}, y0, theta0]

2·nY·np vars: [Xc, L] interleaved, per layer per component (F-order reshape → Xk[k,j])
np vars: h[j] — normalised total height per printed component (0→1, cutoff = Ly·h)
np vars: Mc[j] — membership density per printed component
2 vars: y0, theta0 — global printing-plane y-offset and rotation angle

class ggp.projection.alm_2d.ALM2DMapper(num_components: int, num_layers: int, comp_per_layer: int, layer_height: float, r_gp: float = 0.5, method: str = 'MNA', Ngp: int = 2, **kwargs)[source]

Bases: ProjectionMapper

ALM 2-D projection mapper using the MNA continuous formulation.

All design variables follow the interleaved layout described in the module docstring. The mapper always operates in continuous (MNA) mode — the old per-component 3-var mode is removed.

default_bounds(domain_extents: Tuple[float, ...], num_components: int) → Tuple[numpy.ndarray, numpy.ndarray][source]

Return (lb, ub) arrays for the full design-variable vector.

Parameters

domain_extentstuple of float: (Lx, Ly) for 2-D, (Lx, Ly, Lz) for 3-D.
num_componentsint: Number of geometric primitives.

Returns

lb, ubndarray, ndarray: Lower and upper bounds, each of shape (n_vars,).

forward(x_vars: numpy.ndarray, eval_coords: numpy.ndarray, power_E: float = 1.0, power_V: float = 1.0) → Tuple[numpy.ndarray, numpy.ndarray][source]

Compute density fields.

Parameters

x_varsndarray, shape (n_vars,): Flat design-variable vector (already unscaled).
eval_coordsndarray, shape (n_eval, dim): Element centroid (or evaluation-point) coordinates.
power_Efloat: Penalisation exponent for stiffness density.
power_Vfloat: Penalisation exponent for volume density.

Returns

rho_Endarray, shape (n_elements,): Stiffness density per element.
rho_Vndarray, shape (n_elements,): Volume density per element.

jacobian(x_vars: numpy.ndarray, eval_coords: numpy.ndarray, power_E: float = 1.0, power_V: float = 1.0) → Dict[str, numpy.ndarray][source]

Compute analytic Jacobians.

Returns

dict: {"rho_E": drho_E_dx, "rho_V": drho_V_dx} where each value has shape (n_elements, n_vars).

num_vars() → int[source]: Total number of design variables.

num_vars_per_component() → int[source]: Legacy; returns total vars (not per-component — layout is not uniform).

Optimization Pipeline

Optimization pipeline orchestration.

class ggp.optimization.pipeline.GGPPipeline(spec: ProblemSpec, x0=None, overrides=None, deflation=None)[source]

Bases: object

Orchestrates the entire GGP optimization process.

run() → OptimisationResult[source]

Results containers for GGP optimization.

class ggp.optimization.results.OptimisationResult(problem_name: str, algorithm: str, status: str, iterations: int, objective_value: float, max_constraint_violation: float, design_variables: numpy.ndarray, history: Dict[str, List[float]], objective_name: str = 'compliance', mesh_path: str | None = None, execution_time_s: float | None = None, density_field: numpy.ndarray | None = None, eval_coords: numpy.ndarray | None = None, dim: int = 2)[source]

Bases: object

Container for the results of a GGP optimization run.

algorithm: str

density_field: numpy.ndarray | None = None

design_variables: numpy.ndarray

dim: int = 2

eval_coords: numpy.ndarray | None = None

execution_time_s: float | None = None

history: Dict[str, List[float]]

iterations: int

max_constraint_violation: float

mesh_path: str | None = None

objective_name: str = 'compliance'

objective_value: float

problem_name: str

status: str

to_dict() → Dict[str, Any][source]: Convert the result to a serializable dictionary.

GEMSEO Wrappers

Refactored geometry discipline using projection mappers.

class ggp.gemseo_wrappers.geometry_discipline.GGPGeometryDiscipline(*args: Any, **kwargs: Any)[source]

Bases: Discipline

GEMSEO discipline for GGP geometry projection.

Delegates the core maths to a ProjectionMapper.

Refactored physics discipline using pre-assembled vectors.

class ggp.gemseo_wrappers.physics_discipline.GGPPhysicsDiscipline(*args: Any, **kwargs: Any)[source]

Bases: Discipline

High-performance Physics Discipline using Fast Assembly (petsc4py/scipy/AMJax). Uses directly assembled load vectors and BCs from the FEMDiscretiser.

adjoint_solve(rhs: numpy.ndarray) → numpy.ndarray[source]

Solve the adjoint system K λ = rhs with the same BC-eliminated K.

K is symmetric, so the cached primal factorization (direct solver) is reused directly; iterative/amjax backends re-solve with the same operator. The RHS is zeroed on fixed DOFs to stay consistent with the BC elimination (λ = 0 there).

dstiffness_sensitivity(lam: numpy.ndarray) → numpy.ndarray[source]

Return the element-wise -λᵀ (∂K/∂ρ_e) u term for an adjoint vector lam.

For SIMP K = Σ_e E(ρ_e) kₑ so ∂K/∂ρ_e = dE_drho_e · ke_ref acting on element e’s DOFs only; the implicit-u part of any constraint sensitivity is then dg/dρ_e |_implicit = -λ_eᵀ (dE_drho_e ke_ref) u_e.

Physics Solvers

AMG-preconditioned iterative solver for the FEM stiffness system.

Implements the core approach of the AMJax library (https://github.com/vboussange/AMJax): a PyAMG algebraic multigrid (AMG) hierarchy is built from the assembled stiffness matrix and used as a preconditioner for a conjugate-gradient (CG) solve.

We use scipy’s CG with the PyAMG preconditioner, which avoids the Python ≥ 3.11 requirement of the amjax package while delivering the same convergence behaviour. When JAX is available (and jax_enable_x64 is set), the assembly of the RHS and solution conversion is done via JAX arrays for consistency with future GPU backends.

ggp.physics.amjax_solver.solve_amjax(K: scipy.sparse.spmatrix, f: numpy.ndarray, *, tol: float = 1e-08, maxiter: int = 200, cycle: str = 'V') → numpy.ndarray[source]

Solve K u = f using a PyAMG AMG preconditioner with scipy CG.

Builds a smoothed-aggregation multigrid hierarchy via PyAMG and uses it as a preconditioner for scipy.sparse.linalg.cg(). This reproduces the key algorithmic idea of AMJax: AMG-preconditioned Krylov iteration for FEM linear systems.

Parameters

K:: Assembled, BC-applied global stiffness matrix (any scipy sparse format).
f:: Load vector, shape (n_dofs,).
tol:: Relative residual tolerance for the CG solver.
maxiter:: Maximum CG iterations.
cycle:: AMG cycle type used in the preconditioner — "V", "W", or "F".

Returns

np.ndarray: Displacement vector, shape (n_dofs,), as a NumPy float64 array.