Solver pipeline

This document describes the numerical engine of GenSec: how it goes from a section definition to resistance domains and utilization ratios. Each step is explained with both the engineering meaning and the software mechanism.

Overview: three layers

The solver pipeline has three layers, each building on the previous:

        flowchart LR
    A["Layer 1<br/>FiberSolver<br/><i>strain plane → forces</i>"]
    B["Layer 2<br/>NMDiagram<br/><i>scan configurations → domain</i>"]
    C["Layer 3<br/>VerificationEngine<br/><i>domain + demands → η</i>"]
    A --> B --> C

  • Layer 1 (FiberSolver): the direct problem — given a strain plane, compute the resulting (N, Mx, My). Also solves the inverse problem — given target forces, find the strain plane.

  • Layer 2 (NMDiagram): systematically calls Layer 1 with many strain configurations to build the resistance domain.

  • Layer 3 (VerificationEngine): compares demand points against the domain boundary and computes utilization ratios.

Layer 1: FiberSolver

The FiberSolver lives in solver/integrator.py. It is the computational core of GenSec — everything else is built on top of it.

Direct problem

Given three parameters \((\varepsilon_0, \chi_x, \chi_y)\), the solver computes the internal forces \((N, M_x, M_y)\) in three steps:

        flowchart TD
    P["Input: ε₀, χₓ, χᵧ"]
    S1["Step 1: Compute strain at each fiber<br/>εᵢ = ε₀ + χₓ(yᵢ − y_ref) − χᵧ(xᵢ − x_ref)"]
    S2["Step 2: Compute stress at each fiber<br/>σᵢ = material.stress(εᵢ)"]
    S3["Step 3: Sum forces and moments<br/>N = Σ σᵢ·Aᵢ<br/>Mx = Σ σᵢ·Aᵢ·(yᵢ − y_ref)<br/>My = −Σ σᵢ·Aᵢ·(xᵢ − x_ref)"]
    R["Output: N, Mx, My"]
    P --> S1 --> S2 --> S3 --> R

Step 2 is where the material constitutive law enters: each fiber’s strain is converted to stress using its assigned material’s σ-ε relationship. Different fibers can have different materials (concrete in the bulk, steel at rebar locations, CFRP at an external strip, etc.).

For embedded rebars (bars inside the concrete), the solver subtracts the displaced concrete to avoid double-counting:

\[F_{\text{net},j} = \bigl[\sigma_{\text{rebar}}(\varepsilon_j) - \sigma_{\text{concrete}}(\varepsilon_j)\bigr] \cdot A_{s,j}\]

The subtraction is automatic for all rebars with embedded: true (the default).

Inverse problem (Newton-Raphson)

The inverse problem is: given target forces \((N^*, M_x^*, M_y^*)\), find the strain plane \((\varepsilon_0, \chi_x, \chi_y)\) that produces them.

This is a system of 3 nonlinear equations in 3 unknowns. GenSec solves it with Newton-Raphson iteration with backtracking line search:

        flowchart TD
    A["Initial guess: (ε₀, χₓ, χᵧ)"]
    B["Evaluate integrate_with_tangent<br/>→ (N, Mx, My) + 3×3 tangent K"]
    C{"Residual r = (N−N*, Mx−Mx*, My−My*)<br/>|r| < tolerance?"}
    E["Solve K · Δx = −r"]
    F["Backtracking line search:<br/>find α such that |r(x + α·Δx)| < |r(x)|"]
    G["Update: x ← x + α·Δx"]
    H["Converged!<br/>Return (ε₀, χₓ, χᵧ, N, Mx, My)"]

    A --> B --> C
    C -- No --> E --> F --> G --> B
    C -- Yes --> H

The Jacobian is computed analytically from the tangent stiffness matrix:

\[\mathbf{K} = \sum_i E_{t,i} \, A_i \; \boldsymbol{\varphi}_i \, \boldsymbol{\varphi}_i^T\]

where \(E_{t,i} = d\sigma_i / d\varepsilon_i\) is the tangent modulus at the current strain and \(\boldsymbol{\varphi}_i = [1,\; (y_i - y_{\text{ref}}),\; -(x_i - x_{\text{ref}})]^T\) is the shape-function vector. This is computed in a single pass together with the internal forces by integrate_with_tangent(), halving the cost per iteration compared to the former finite-difference approach.

The line search halves the step size up to 15 times if the full Newton step would increase the residual.

For uniaxial bending (\(\chi_y = 0\)), the system reduces to 2×2 (the upper-left sub-block of \(\mathbf{K}\)), which converges faster. GenSec detects this case automatically.

Layer 2: NMDiagram

The NMDiagram class lives in solver/capacity.py. It uses the FiberSolver to systematically explore the space of possible strain configurations.

N-M diagram (uniaxial)

The 2D resistance domain \((N, M_x)\) is generated by scanning pairs of edge strains \((\varepsilon_{\text{bottom}}, \varepsilon_{\text{top}})\) across the full range of material limits:

        flowchart TD
    A["Collect strain limits:<br/>ε_cu2 (concrete crush),<br/>ε_su (steel ultimate)"]
    B["Branch 1: top edge at ε_cu2<br/>sweep bottom from ε_su to ε_cu2"]
    C["Branch 2: bottom edge at ε_cu2<br/>sweep top from ε_su to ε_cu2"]
    D["Branch 3: uniform tension<br/>(both edges positive)"]
    E["Branch 4: near-uniform compression<br/>(both edges near ε_cu2)"]
    F["For each (ε_bot, ε_top):<br/>convert to (ε₀, χₓ)<br/>evaluate direct problem<br/>collect (N, Mx)"]
    G["Result: point cloud<br/>on the N-M boundary"]

    A --> B & C & D & E
    B & C & D & E --> F --> G

This is not the classic EC2 pivot method (which identifies specific pivot points A, B, C). GenSec uses a parametric scan that covers the same configuration space but is more general — it works with any material, not just EC2 parabola-rectangle.

The scan covers both curvature directions (positive and negative \(\chi_x\)), which is essential for producing a symmetric diagram.

3D resistance surface (biaxial)

For biaxial bending, the curvature direction is parameterized by an angle \(\theta\):

\[\chi_x = \chi \cos\theta, \qquad \chi_y = \chi \sin\theta\]

The scan loops over \(\theta \in [0, 2\pi)\) in n_angles steps, and for each angle performs the same edge-strain scan as in the uniaxial case. This produces a 3D point cloud \((N, M_x, M_y)\).

Mx-My contour at fixed N

A horizontal slice of the 3D surface. For each curvature direction \(\theta\), GenSec scans strain configurations, collects all \((N, M_x, M_y)\) points, and interpolates at \(N = N_{\text{fixed}}\) to find the moment point on the contour.

Among all crossing points at the target N, the one with maximum moment magnitude is kept. The result is a closed curve in the \((M_x, M_y)\) plane.

Moment-curvature (M-χ)

At a fixed axial force \(N\), GenSec scans curvature \(\chi\) from 0 to ultimate (and from 0 to −ultimate for the other direction). At each curvature step:

  1. Solve for \(\varepsilon_0\) such that \(N(\varepsilon_0, \chi) = N_{\text{fixed}}\) (1D Newton).

  2. Compute the moment \(M(\varepsilon_0, \chi)\).

  3. Check for key events: cracking (concrete tensile strain exceeds \(\varepsilon_{cr} = f_{ctm}/E_{cm}\)), first yield (any rebar strain exceeds \(\varepsilon_{yd}\)), ultimate (any material reaches its strain limit).

Layer 3: VerificationEngine

The VerificationEngine lives in solver/check.py. It wraps the resistance domain and provides four utilization ratio methods.

Architecture

        classDiagram
    class VerificationEngine {
        -domain: DomainChecker
        -nm_gen: NMDiagram
        -_contour_cache: dict
        +check_demand(demand) dict
        +check_combination(combo, demand_db) dict
        +check_envelope(envelope, demand_db) dict
    }
    class DomainChecker {
        -hull: ConvexHull
        -N_range: float
        +eta_3D(N, Mx, My) float
        +eta_path(base, target) float
        +is_inside(N, Mx, My) bool
    }
    class MxMyContour {
        -hull: ConvexHull (2D)
        +eta_2D(Mx, My) float
        +eta_path_2D(base, target) float
    }

    VerificationEngine --> DomainChecker : 3D operations
    VerificationEngine --> MxMyContour : 2D operations (cached)

  • DomainChecker: wraps the 3D ConvexHull. Handles eta_3D (ray from origin) and eta_path (ray from arbitrary base). Always available — the 3D hull is generated for any biaxial section.

  • MxMyContour: wraps a 2D ConvexHull of the Mx-My contour at a specific N level. Handles eta_2D and eta_path_2D. Generated on demand and cached by the engine.

  • VerificationEngine: the orchestrator. Reads the output flags, resolves demand references, accumulates combination stages, and dispatches to the appropriate checker.

The ray-casting algorithm

All four η types use the same geometric primitive:

  1. Define a ray from base \(\mathbf{B}\) through target \(\mathbf{T}\).

  2. Find the first intersection \(\mathbf{R}\) of the ray with the ConvexHull boundary (by solving a linear system against each hull facet).

  3. Compute \(\eta = |\mathbf{T} - \mathbf{B}| \;/\; |\mathbf{R} - \mathbf{B}|\).

        flowchart LR
    B["Base B"] -->|"ray direction"| T["Target T"]
    T -.->|"continues to"| R["Boundary R"]
    style B fill:#4CAF50,color:#fff
    style T fill:#F44336,color:#fff
    style R fill:#2196F3,color:#fff

If \(\eta < 1\), the target is closer to the base than the boundary — the section has reserve capacity. If \(\eta > 1\), the target is beyond the boundary — the section fails.

Staged combinations

For staged loading (e.g., gravity → seismic), the engine accumulates forces stage by stage:

        flowchart TD
    O["Origin (0,0,0)"]
    S0["Stage 0: Gravity<br/>S₀ = 1.0·G + 0.3·Q"]
    S1["Stage 1: Seismic<br/>S₁ = S₀ + 1.0·Ex"]
    R0["η_3D of S₀<br/>(ray from O to S₀)"]
    R1["η_path of S₁<br/>(ray from S₀ to S₁)"]

    O -->|"η_3D"| S0
    S0 -->|"η_path"| S1
    S0 --> R0
    S1 --> R1

The engineering meaning of \(\eta_{\text{path}}\) for the seismic stage is: “the gravity load is certain; how much of the seismic capacity is being used?” This is the relevant question for seismic design, because the gravity load does not vary — only the seismic action is the “variable” that might exhaust capacity.

Material constitutive laws

Each material class inherits from the abstract Material base and implements four methods:

  • stress(eps) — scalar: one strain in, one stress out.

  • stress_array(eps) — vectorized: array in, array out (fast path). Accepts arrays of any shape (1-D, 2-D, …).

  • tangent(eps) — scalar tangent modulus \(E_t = d\sigma/d\varepsilon\).

  • tangent_array(eps) — vectorized tangent modulus (any shape).

        flowchart LR
    E["ε (strain)"] --> M["Material.stress(ε)"] --> S["σ (stress)"]
    E --> T["Material.tangent(ε)"] --> Et["Eₜ (tangent)"]

The fiber solver calls stress_array and tangent_array on entire fiber arrays at once, leveraging NumPy vectorization for speed. When numba is installed, the built-in Concrete and Steel classes use JIT-compiled kernels for these methods, providing a further ~2–3× speed-up.

Concrete (parabola-rectangle, EC2 3.1.7):

\[\begin{split}\sigma_c(\varepsilon) = \begin{cases} 0 & \varepsilon > 0 \quad \text{(no tension)} \\ -f_{cd}\!\left[1 - \left(1 - \frac{\varepsilon} {\varepsilon_{c2}}\right)^{n}\right] & \varepsilon_{c2} \le \varepsilon \le 0 \\ -f_{cd} & \varepsilon_{cu2} \le \varepsilon < \varepsilon_{c2} \\ 0 & \varepsilon < \varepsilon_{cu2} \quad \text{(crushed)} \end{cases}\end{split}\]

Steel (elastic-plastic with optional hardening):

\[\begin{split}\sigma_s(\varepsilon) = \begin{cases} E_s \, \varepsilon & |\varepsilon| \le \varepsilon_{yd} \\ f_{yd} + (k-1)\,f_{yd}\, \frac{|\varepsilon| - \varepsilon_{yd}} {\varepsilon_{su} - \varepsilon_{yd}} & \varepsilon_{yd} < |\varepsilon| \le \varepsilon_{su} \\ 0 & |\varepsilon| > \varepsilon_{su} \end{cases}\end{split}\]

TabulatedMaterial: piecewise-linear interpolation of user-supplied (ε, σ) data points. Supports any shape, including descending branches (Kent-Park concrete, softening, etc.). The tangent modulus is the slope of each linear segment.

Performance architecture

GenSec uses three complementary strategies to achieve high computational throughput on the fiber integration pipeline.

Batch integration

The capacity generators (N-M diagram, 3-D surface, Mx-My contour) evaluate tens of thousands of strain configurations. Rather than calling integrate() in a Python loop, they use integrate_batch(), which operates on 1-D arrays of \((\varepsilon_0, \chi_x, \chi_y)\) and produces 1-D arrays of \((N, M_x, M_y)\).

Internally, the strain field is built as a 2-D array of shape \((n_{\text{configs}}, n_{\text{fibers}})\) and passed through the constitutive law in one NumPy call. This eliminates the per-call Python overhead and enables SIMD vectorization at the C level.

Mega-batch with chunking

For biaxial generators (generate_biaxial, generate_mx_my), configurations from all curvature directions are concatenated into a single flat array before integration. This avoids the overhead of calling integrate_batch once per angle (e.g. 144 calls of 550 configs each) and replaces it with a few large chunked calls (e.g. 2 calls of 40 000 configs each).

The chunk size is bounded to limit peak memory usage:

\[n_{\text{chunk}} = \max\!\left(2000,\; \left\lfloor\frac{50 \times 10^6}{n_{\text{fibers}}} \right\rfloor\right)\]

For a section with 5 000 fibers, this gives chunks of 10 000 configs (~400 MB peak), reducing the number of batch calls from 144 to ~8.

Optional Numba JIT

When numba is installed, the element-wise stress and tangent kernels of Concrete and Steel are JIT-compiled to native code. This provides a ~2–3× speed-up on large fiber arrays where NumPy’s Python-to-C dispatch overhead becomes significant.

Numba is optional: if not installed, the pure-NumPy fallback is used transparently. Install with pip install gensec[fast] or uv sync --all-extras.

Section meshing

The GenericSection class (geometry/geometry.py) takes a Shapely polygon and produces a fiber mesh:

        flowchart TD
    P["Shapely Polygon<br/>(exterior + holes)"]
    M{"Mesh method?"}
    G["Grid meshing:<br/>rectangular grid cells<br/>clipped to polygon boundary"]
    T["Triangle meshing:<br/>constrained Delaunay<br/>triangulation"]
    F["Fiber arrays:<br/>x_fibers, y_fibers, A_fibers"]
    R["Rebar arrays:<br/>x_rebars, y_rebars, A_rebars"]

    P --> M
    M -- grid --> G --> F
    M -- triangle --> T --> F
    P --> R

Each grid cell is intersected with the polygon. Partial cells at the boundary receive the area of the intersection, not the full cell area. This gives a good area representation even with coarse meshes.

When explicit grid dimensions are provided via n_grid_x / n_grid_y (used internally by the RectSection factory function), they override the mesh_size-based computation. This allows anisotropic grids when needed.

The RectSection convenience function creates a rectangular polygon and calls GenericSection internally. It derives an isotropic mesh size from n_fibers_y (i.e., \(s = H / n_y\)), ensuring approximately square cells. When n_fibers_x > 1 is specified explicitly, it is passed as n_grid_x to override the x-resolution.