Function bodies 335 total

class Colors:
    RED = "\033[91m"
    YELLOW = "\033[93m"
    GREEN = "\033[92m"
    BLUE = "\033[94m"
    CYAN = "\033[96m"
    BOLD = "\033[1m"
    DIM = "\033[2m"
    RESET = "\033[0m"

def risk_badge(rating: str) -> str:
    if rating == "RED":
        return f"{Colors.RED}🔴 RED{Colors.RESET}"
    elif rating == "YELLOW":
        return f"{Colors.YELLOW}🟡 YELLOW{Colors.RESET}"
    else:
        return f"{Colors.GREEN}🟢 GREEN{Colors.RESET}"

def parse_description(text: str) -> dict:
    """Extract bridge parameters from natural language description.
    
    This is deliberately simple — it extracts what it can from patterns
    and leaves the rest for the user to confirm or the tools to infer.
    """
    import re
    
    params = {
        "description": text.strip(),
        "bridge_type": None,
        "spans": [],
        "num_girders": None,
        "girder_spacing_ft": None,
        "deck_width_ft": None,
        "girder_depth_in": None,
        "erection_method": None,
        "location": {},
        "water_crossing": False,
        "constraints": [],
    }
    
    lower = text.lower()
    
    # Bridge type
    type_patterns = {
        "steel_plate_girder": r"steel\s+(?:plate\s+)?girder",
        "prestressed_i": r"prestress(?:ed)?\s+(?:concrete\s+)?i[- ]?girder",
        "concrete_box": r"(?:cip|cast.in.place|concrete)\s+box\s+girder",
        "segmental_box": r"(?:segmental|pt\s+segmental)\s+box",
        "s

def _print_force_summary(analysis: dict, superstructure: dict) -> None:
    """Print a compact shear/moment summary table after analysis."""
    envelopes = analysis.get("envelopes", {})
    ml_envelopes = analysis.get("moving_load_envelopes", {})
    reactions = analysis.get("reactions", {})
    displacements = analysis.get("displacements", {})
    
    if not envelopes:
        return
    
    # Gather gravity envelope peaks
    max_pos_M = 0.0  # max positive moment (sagging)
    max_neg_M = 0.0  # max negative moment (hogging)
    max_V = 0.0
    ctrl_pos_elem = ""
    ctrl_neg_elem = ""
    ctrl_V_elem = ""
    
    for etag, env in envelopes.items():
        mz = env.get("Mz_i", {})
        vy = env.get("Vy_i", {})
        if isinstance(mz, dict):
            mx = mz.get("max", 0)
            mn = mz.get("min", 0)
            if mx > max_pos_M:
                max_pos_M = mx; ctrl_pos_elem = etag
            if abs(mn) > abs(max_neg_M):
                max_neg_M = mn; ctrl_neg_el

def run_pipeline(params: dict, output_dir: Path | None = None, verbose: bool = False, dump_script: bool = False):
    """Run the full NLB pipeline from parsed parameters."""
    
    t0 = time.time()
    results = {}
    
    # ── Step 1: Site Recon ────────────────────────────────────────────
    status("🔍", "Site reconnaissance...")
    
    lat = params.get("location", {}).get("lat")
    lon = params.get("location", {}).get("lon")
    
    if lat and lon:
        try:
            from nlb.tools.site_recon import run_site_recon
            site = run_site_recon(lat, lon, params["description"])
            site_dict = site.to_dict() if hasattr(site, "to_dict") else site
            results["site"] = site_dict
            
            sdc = site_dict.get("seismic", {}).get("sdc", "?")
            wind = site_dict.get("wind", {}).get("v_ult", "?")
            scour = "scour zone" if site_dict.get("scour", {}).get("water_crossing") else "no scour"
            state = params["location"].g

def _build_superstructure_params(params: dict) -> dict:
    """Build superstructure tool parameters from parsed description."""
    spans = params.get("spans", [100.0])
    btype = params.get("bridge_type", "steel_plate_girder")
    
    type_map = {
        "steel_plate_girder": "steel_plate_girder_composite",
        "prestressed_i": "prestressed_i_girder",
        "concrete_box": "cip_box_girder",
        "segmental_box": "segmental_box_girder",
        "steel_truss": "steel_truss",
        "concrete_slab": "concrete_slab",
        "arch": "arch",
        "cable_stayed": "cable_stayed",
    }
    
    # Use single-line model (1 girder) for analysis.
    # Multi-girder models cause force amplification from penalty constraints.
    # Distribution factors handle transverse load distribution.
    p = {
        "bridge_type": type_map.get(btype, "steel_plate_girder_composite"),
        "span_lengths_ft": spans,
        "num_girders": 1,  # Single-line model; DFs applied in post-processin

def _build_foundation_params(params: dict, index: int, total: int, site: dict) -> dict:
    """Build foundation parameters for a given support."""
    is_abutment = (index == 0 or index == total - 1)
    
    return {
        "foundation_type": "spread_footing" if is_abutment else "drilled_shaft",
        "params": {
            "diameter_ft": 4.0 if is_abutment else 7.0,
            "length_ft": 15.0 if is_abutment else 60.0,
            "width_ft": 20.0 if is_abutment else None,
            "depth_ft": 6.0 if is_abutment else None,
        },
        "site_profile": site,
    }

def _build_substructure_params(params: dict, index: int, total: int) -> dict:
    """Build substructure parameters for a given support."""
    is_abutment = (index == 0 or index == total - 1)
    
    if is_abutment:
        return {
            "sub_type": "seat_abutment",
            "seat_width_ft": 4.0,
            "backwall_height_ft": 7.0,
            "num_bearings": params.get("num_girders", 5),
            "bearing_spacing_ft": params.get("girder_spacing_ft", 8.0),
        }
    else:
        return {
            "sub_type": "multi_column_bent",
            "num_columns": 2,
            "column_diameter_ft": 7.0,
            "column_height_ft": 25.0,
            "column_spacing_ft": 20.0,
            "cap_width_ft": 4.0,
            "cap_depth_ft": 5.0,
            "fc_ksi": 4.0,
            "fy_ksi": 60.0,
        }

def _build_bearing_params(params: dict, index: int, total: int) -> dict:
    """Build bearing parameters for a given support."""
    is_abutment = (index == 0 or index == total - 1)
    
    return {
        "bearing_type": "elastomeric" if is_abutment else "ptfe_sliding",
        "vertical_capacity_kips": 500.0,
    }

def _build_load_params(params: dict, site: dict) -> dict:
    """Build load generation parameters."""
    from nlb.tools.loads import BridgeGeometry
    
    spans = params.get("spans", [100.0])
    n_girders = params.get("num_girders", 5)
    spacing = params.get("girder_spacing_ft", 8.0)
    deck_w = params.get("deck_width_ft") or (n_girders * spacing)
    depth = params.get("girder_depth_in") or 48.0
    
    # BridgeGeometry.span_ft expects a single float (max span for load generation)
    max_span = max(spans) if spans else 100.0
    
    geom = BridgeGeometry(
        span_ft=max_span,
        girder_spacing_ft=spacing,
        num_girders=n_girders,
        deck_width_ft=deck_w,
        girder_depth_in=depth,
    )
    
    return {
        "geom": geom,
        "site_profile": site,
    }

def _model_to_dict(obj) -> dict:
    """Convert a dataclass or object to dict."""
    if isinstance(obj, dict):
        return obj
    if hasattr(obj, "__dataclass_fields__"):
        import dataclasses
        return dataclasses.asdict(obj)
    if hasattr(obj, "to_dict"):
        return obj.to_dict()
    if hasattr(obj, "__dict__"):
        return {k: v for k, v in obj.__dict__.items() if not k.startswith("_")}
    return {}

def _default_site(params: dict) -> dict:
    """Conservative default site profile when recon fails."""
    return {
        "coordinates": params.get("location", {"lat": 40.0, "lon": -90.0}),
        "seismic": {"pga": 0.10, "ss": 0.25, "s1": 0.10, "sds": 0.20, "sd1": 0.13, "site_class": "D", "sdc": "B"},
        "wind": {"v_ult": 115, "exposure": "C"},
        "thermal": {"t_min": -10, "t_max": 110, "delta_t": 120},
        "scour": {"water_crossing": params.get("water_crossing", False), "design_flood": "Q100", "check_flood": "Q500"},
        "frost_depth_ft": 3.0,
        "soil": {"site_class": "D", "description": "Stiff soil (default)"},
        "climate_zone": "cold",
        "layers": [
            {"soil_type": "stiff_clay", "top_depth_ft": 0, "thickness_ft": 15, "su_ksf": 1.5, "gamma_pcf": 120},
            {"soil_type": "sand", "top_depth_ft": 15, "thickness_ft": 25, "phi_deg": 35, "gamma_pcf": 125, "N_spt": 30},
            {"soil_type": "stiff_clay", "top_depth_ft": 40, "thic

def main():
    parser = argparse.ArgumentParser(
        prog="nlb",
        description="Red Team Your Bridge — Natural Language → Nonlinear FEA → Adversarial Analysis",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  nlb "3-span steel plate girder, 315-420-315 ft, Illinois"
  nlb --input bridge.txt --output ./report/
  nlb site-recon --lat 42.28 --lon -89.09
        """,
    )
    
    parser.add_argument("description", nargs="*", help="Bridge description in natural language")
    parser.add_argument("--input", "-i", help="Read description from file")
    parser.add_argument("--output", "-o", default="./nlb-output", help="Output directory (default: ./nlb-output)")
    parser.add_argument("--verbose", "-v", action="store_true", help="Show detailed output and tracebacks")
    parser.add_argument("--json", action="store_true", help="Output raw JSON instead of pretty print")
    parser.add_argument("--parse-only", action="store_true", help="O

def convergence_handler(tol: float = 1.0e-8, max_iter: int = 100) -> bool:
    """Set up Newton-based convergence fallback chain.

    Attempts progressively more robust solution algorithms:
    1. Newton-Raphson (fastest convergence near solution)
    2. Modified Newton (cheaper iterations, slower convergence)
    3. BFGS (quasi-Newton, good for ill-conditioned problems)
    4. Broyden (quasi-Newton variant, last resort)

    Call this BEFORE running analysis steps. If Newton fails during analysis,
    call try_algorithms() to attempt fallback solutions.

    Args:
        tol:      Convergence tolerance (norm of unbalanced force).
                  Default: 1e-8. Use 1e-6 for pushover, 1e-10 for modal.
        max_iter: Maximum iterations per step. Default: 100.

    Returns:
        True if test and algorithm were set successfully.

    Reference:
        OpenSees analysis command documentation.
        Scott, M.H. (2011). "Numerical Integration Options for the
        Force-Based B

def try_algorithms(tol: float = 1.0e-8, max_iter: int = 200) -> int:
    """Attempt analysis step with fallback algorithm chain.

    Tries Newton → ModifiedNewton → BFGS → Broyden until convergence
    is achieved or all fail.

    Args:
        tol:      Convergence tolerance.
        max_iter: Maximum iterations for fallback algorithms.

    Returns:
        0 if successful, negative if all algorithms fail.

    Usage:
        >>> ok = ops.analyze(1)
        >>> if ok != 0:
        ...     ok = try_algorithms()
    """
    # Try progressively more robust algorithms, tests, and tolerances
    configs = [
        ('NormUnbalance', 'Newton', [], tol, max_iter),
        ('NormUnbalance', 'ModifiedNewton', [], tol, max_iter),
        ('NormUnbalance', 'Newton', [], tol * 100, max_iter),
        ('EnergyIncr', 'Newton', [], tol, max_iter),
        ('EnergyIncr', 'ModifiedNewton', [], tol, max_iter),
        ('NormUnbalance', 'BFGS', [], tol * 100, max_iter),
        ('NormUnbalance', 'Bro

def gravity_analysis(steps: int = 10, tol: float = 1.0e-8) -> int:
    """Run load-controlled static gravity analysis.

    Applies all currently defined loads in `steps` equal increments.
    Uses load-control integrator with fallback convergence handling.

    Typical usage:
        1. Define model (nodes, elements, materials)
        2. Apply gravity loads (nodal, elemental, self-weight)
        3. Call gravity_analysis()
        4. Set loads constant: ops.loadConst('-time', 0.0)

    Args:
        steps: Number of load increments. Default: 10.
               More steps = better convergence for nonlinear problems.
        tol:   Convergence tolerance. Default: 1e-8.

    Returns:
        0 if successful, negative if analysis fails.

    Reference:
        AASHTO LRFD 3.4.1: Load Factors and Load Combinations.
    """
    ops.system('BandGeneral')
    ops.numberer('RCM')
    ops.constraints('Plain')
    ops.integrator('LoadControl', 1.0 / steps)
    convergence_handler(tol)
    ops.a

class PushoverResult:
    """Results from pushover analysis.

    Attributes:
        displacements: List of displacement values at control node (inches).
        base_shear:    List of base shear values (kips).
        steps_completed: Number of steps successfully completed.
        converged:     True if all steps completed.
    """
    displacements: List[float] = field(default_factory=list)
    base_shear: List[float] = field(default_factory=list)
    steps_completed: int = 0
    converged: bool = False

def pushover_analysis(node: int, dof: int, target_disp: float,
                      steps: int = 100, tol: float = 1.0e-6) -> PushoverResult:
    """Run displacement-controlled pushover analysis.

    Monotonically pushes a control node to a target displacement while
    recording the force-displacement response. Essential for:
    - Seismic capacity evaluation
    - Ductility assessment
    - Plastic hinge characterization

    Args:
        node:        Control node tag.
        dof:         DOF to push (1=x, 2=y, 3=z).
        target_disp: Target displacement (inches).
        steps:       Number of displacement increments. Default: 100.
        tol:         Convergence tolerance. Default: 1e-6.

    Returns:
        PushoverResult with force-displacement data.

    Reference:
        AASHTO SGS 4.8: Displacement Capacity Verification.
        Caltrans SDC 5.2.3: Pushover Analysis.
    """
    result = PushoverResult()
    d_incr = target_disp / steps

    ops.system('BandGeneral')

def _adaptive_step(node: int, dof: int, d_incr: float,
                   tol: float, subdivisions: int = 10) -> int:
    """Subdivide a failed pushover step into smaller increments.

    Args:
        node:         Control node.
        dof:          DOF direction.
        d_incr:       Original increment size.
        tol:          Convergence tolerance.
        subdivisions: Number of sub-steps.

    Returns:
        0 if all sub-steps succeed, -1 otherwise.
    """
    sub_incr = d_incr / subdivisions
    for j in range(subdivisions):
        ops.integrator('DisplacementControl', node, dof, sub_incr)
        ok = try_algorithms(tol)
        if ok != 0:
            return -1
    # Restore original increment
    ops.integrator('DisplacementControl', node, dof, d_incr)
    return 0

class ModalResult:
    """Results from modal / response spectrum analysis.

    Attributes:
        periods:          Natural periods (seconds).
        frequencies:      Natural frequencies (Hz).
        mode_shapes:      Dict of {mode: {node: [displacements]}}.
        spectral_accels:  Spectral accelerations at each period (g).
        modal_forces:     Combined forces per CQC.
    """
    periods: List[float] = field(default_factory=list)
    frequencies: List[float] = field(default_factory=list)
    mode_shapes: Dict[int, Dict[int, List[float]]] = field(default_factory=dict)
    spectral_accels: List[float] = field(default_factory=list)
    modal_forces: Dict[str, float] = field(default_factory=dict)

def response_spectrum(damping: float, periods: List[float],
                      accels: List[float],
                      num_modes: int = 10,
                      direction: int = 1) -> ModalResult:
    """Perform modal analysis and response spectrum combination.

    Steps:
    1. Eigenvalue analysis to find natural periods/mode shapes
    2. Read spectral acceleration at each period
    3. Compute modal responses
    4. Combine using CQC (Complete Quadratic Combination)

    Args:
        damping:    Modal damping ratio (e.g., 0.05 for 5%).
        periods:    Design spectrum periods (seconds).
        accels:     Design spectrum accelerations (g) at each period.
        num_modes:  Number of modes to extract. Default: 10.
        direction:  Excitation direction (1=x, 2=y, 3=z). Default: 1.

    Returns:
        ModalResult with periods, mode shapes, and combined responses.

    Reference:
        AASHTO LRFD 4.7.4.3: Multimode Spectral Analysis.
        AASHTO SGS 4.3.3: Respo

def _interpolate_spectrum(T: float, periods: List[float],
                          accels: List[float]) -> float:
    """Linearly interpolate spectral acceleration at period T.

    Args:
        T:       Target period (seconds).
        periods: Spectrum periods (must be sorted ascending).
        accels:  Spectrum accelerations (g).

    Returns:
        Interpolated spectral acceleration (g).
    """
    if not periods or not accels:
        return 0.0
    if T <= periods[0]:
        return accels[0]
    if T >= periods[-1]:
        return accels[-1]

    for i in range(len(periods) - 1):
        if periods[i] <= T <= periods[i + 1]:
            # Linear interpolation
            frac = (T - periods[i]) / (periods[i + 1] - periods[i])
            return accels[i] + frac * (accels[i + 1] - accels[i])

    return accels[-1]

def cqc_combination(responses: List[float], periods: List[float],
                    damping: float) -> float:
    """Combine modal responses using Complete Quadratic Combination (CQC).

    Args:
        responses: Peak modal response for each mode.
        periods:   Natural period for each mode (seconds).
        damping:   Modal damping ratio.

    Returns:
        Combined response.

    Reference:
        Der Kiureghian, A. (1981). "A Response Spectrum Method for
        Random Vibration Analysis of MDF Systems." Earthquake Eng. & Struct. Dyn.
    """
    n = len(responses)
    total = 0.0

    for i in range(n):
        for j in range(n):
            if periods[i] <= 0 or periods[j] <= 0:
                rho = 0.0
            else:
                beta = periods[j] / periods[i]  # frequency ratio
                zeta = damping
                # CQC correlation coefficient
                num = 8.0 * zeta ** 2 * (1.0 + beta) * beta ** 1.5
                den = ((1.0 - beta ** 2)

class TimeHistoryResult:
    """Results from time history analysis.

    Attributes:
        time:              Time vector (seconds).
        steps_completed:   Number of steps completed.
        converged:         True if all steps completed.
        dt_used:           Actual time step used.
    """
    time: List[float] = field(default_factory=list)
    steps_completed: int = 0
    converged: bool = False
    dt_used: float = 0.0

def time_history(dt: float, record: List[float], damping: float = 0.05,
                 tsTag: int = 1, patternTag: int = 1,
                 direction: int = 1, tol: float = 1.0e-8,
                 scale: float = 386.4) -> TimeHistoryResult:
    """Run Newmark transient time history analysis.

    Applies a ground motion acceleration record using the Newmark-beta
    method (average acceleration: gamma=0.5, beta=0.25).

    Args:
        dt:          Time step of the input record (seconds).
        record:      Acceleration time history values (g).
        damping:     Rayleigh damping ratio. Default: 0.05 (5%).
        tsTag:       Time series tag. Default: 1.
        patternTag:  Load pattern tag. Default: 1.
        direction:   Excitation direction (1=x, 2=y, 3=z). Default: 1.
        tol:         Convergence tolerance. Default: 1e-8.
        scale:       Scale factor to convert g to in/s² (386.4 for KIS).
                     Default: 386.4.

    Returns:
        TimeHistoryRes

class ConstructionStage:
    """Definition of a construction stage.

    Attributes:
        name:           Descriptive name (e.g., "Erect Girders", "Pour Deck").
        activate_elements:  Element tags to activate in this stage.
        activate_loads:     Load pattern tags to apply.
        deactivate_elements: Element tags to remove.
        steps:          Number of analysis steps for this stage.
        description:    Engineering description for reporting.
    """
    name: str
    activate_elements: List[int] = field(default_factory=list)
    activate_loads: List[int] = field(default_factory=list)
    deactivate_elements: List[int] = field(default_factory=list)
    steps: int = 10
    description: str = ""

def staged_construction(stages: List[ConstructionStage],
                        tol: float = 1.0e-8) -> Dict[str, int]:
    """Run sequential staged construction analysis.

    Activates/deactivates elements and loads in sequence, running a
    gravity-type analysis at each stage. Critical for:
    - Steel erection sequence
    - Deck pour sequence (wet concrete loads)
    - Barrier/overlay placement
    - PT stressing sequence
    - Shoring removal

    After each stage, loads are set constant (ops.loadConst) so
    subsequent stages see cumulative effects.

    Args:
        stages: Ordered list of ConstructionStage definitions.
        tol:    Convergence tolerance. Default: 1e-8.

    Returns:
        Dict of {stage_name: status_code}. 0 = success, negative = failure.

    Reference:
        AASHTO LRFD 5.12.5: Segmental Construction.
        AASHTO LRFD C3.4.2: Construction Loads.
    """
    results: Dict[str, int] = {}

    for stage in stages:
        # Activate elements (in O

def geometric_transform(tag: int, transform_type: str,
                        *vecxz: float) -> int:
    """Define a geometric transformation for beam-column elements.

    The transformation maps element local coordinates to global coordinates
    and determines how geometric nonlinearity is handled.

    Args:
        tag:             Transformation tag.
        transform_type:  One of:
                         - "Linear"        — small displacement, no P-delta
                         - "PDelta"        — includes P-delta effects
                         - "Corotational"  — large displacement (co-rotational)
        *vecxz:          Components of the vector in the local x-z plane
                         (used to define the local coordinate system).
                         For 2D: not needed (can be omitted).
                         For 3D: typically (0, 0, 1) for vertical elements or
                                 (0, 1, 0) for horizontal elements.

    Returns:
        Transfo

def beam_column(tag: int, nodes: Tuple[int, int], section: int,
                transform: int, integration: str = 'Lobatto',
                np: int = 5) -> int:
    """Create a displacement-based beam-column element.

    Uses dispBeamColumn for displacement-based formulation, which is
    recommended for most bridge analysis applications. The element uses
    distributed plasticity with fiber sections at integration points.

    Args:
        tag:          Element tag.
        nodes:        Tuple of (iNode, jNode) — end node tags.
        section:      Section tag (fiber section from sections module).
        transform:    Geometric transformation tag.
        integration:  Integration scheme. Options:
                      - 'Lobatto'       — Gauss-Lobatto (default, recommended)
                      - 'Legendre'      — Gauss-Legendre
                      - 'Radau'         — Gauss-Radau
                      - 'NewtonCotes'   — Newton-Cotes
                      - 'HingeRadau'    

def truss_element(tag: int, nodes: Tuple[int, int], area: float,
                  material: int) -> int:
    """Create a co-rotational truss element.

    Uses corotTruss for large-displacement truss analysis, suitable for
    cable elements, bracing, and truss bridges. The co-rotational
    formulation accounts for geometric nonlinearity.

    Args:
        tag:      Element tag.
        nodes:    Tuple of (iNode, jNode).
        area:     Cross-sectional area (in²).
        material: Uniaxial material tag.

    Returns:
        Element tag.

    Reference:
        AASHTO LRFD 4.6.2.7: Truss Bridges.
    """
    ops.element('corotTruss', tag, *nodes, area, material)
    return tag

def zero_length(tag: int, nodes: Tuple[int, int],
                materials: List[int], directions: List[int]) -> int:
    """Create a zero-length element for bearings, springs, or releases.

    Zero-length elements connect two nodes at the same location with
    specified material behavior in each DOF. Used extensively for:
    - Elastomeric bearing pads (shear + axial)
    - Soil springs (p-y, t-z, q-z)
    - Moment releases (pins)
    - Friction elements

    Args:
        tag:         Element tag.
        nodes:       Tuple of (iNode, jNode) — must be at same coordinates.
        materials:   List of material tags, one per direction.
        directions:  List of DOF directions (1=x, 2=y, 3=z, 4=rx, 5=ry, 6=rz).
                     Must match length of materials.

    Returns:
        Element tag.

    Raises:
        ValueError: If materials and directions have different lengths.

    Reference:
        AASHTO LRFD 14.7: Bearings.
        OpenSees zeroLength command documentation

def shell_element(tag: int, nodes: Tuple[int, int, int, int],
                  thickness: float, material: int) -> int:
    """Create a ShellMITC4 element for deck or wall modeling.

    MITC4 (Mixed Interpolation of Tensorial Components) formulation avoids
    shear locking and is suitable for both thin and moderately thick shells.

    Args:
        tag:        Element tag.
        nodes:      Tuple of 4 node tags (counter-clockwise ordering).
        thickness:  Shell thickness (inches).
        material:   nDMaterial tag (e.g., ElasticIsotropic or PlateFromPlaneStress).

    Returns:
        Element tag.

    Note:
        The material must be an nDMaterial (2D), not a uniaxialMaterial.
        For bridge decks, typically use:
            ops.nDMaterial('ElasticIsotropic', matTag, E, nu)
            ops.section('PlateFiber', secTag, matTag, thickness)
        Then pass secTag to this function.

    Reference:
        AASHTO LRFD 4.6.3.3: Refined Methods (FEM).
        Dvorkin, E.N

class ConcreteProperties:
    """Computed concrete material properties.

    Attributes:
        fc:    Compressive strength (ksi, positive value)
        Ec:    Elastic modulus (ksi)
        fr:    Modulus of rupture (ksi)
        eps_0: Strain at peak stress (in/in, positive value)
        eps_cu: Ultimate crushing strain (in/in, positive value)
    """
    fc: float
    Ec: float
    fr: float
    eps_0: float
    eps_cu: float

def concrete_defaults(fc_ksi: float) -> ConcreteProperties:
    """Auto-compute concrete properties from f'c.

    Per AASHTO LRFD 5.4.2.4 and ACI 318-19 Table 19.2.2.1.

    Args:
        fc_ksi: Specified compressive strength f'c in ksi (positive).

    Returns:
        ConcreteProperties with Ec, fr, eps_0, eps_cu.

    Example:
        >>> props = concrete_defaults(4.0)  # 4 ksi concrete
        >>> round(props.Ec, 0)  # ~3644 ksi
        3644.0
    """
    fc = abs(fc_ksi)
    # Unit weight of normal-weight concrete: 0.150 kcf = 150 pcf
    wc = 0.150  # kcf

    # AASHTO LRFD Eq. 5.4.2.4-1: Ec = 33000 * K1 * wc^1.5 * sqrt(f'c)
    # K1 = 1.0 (correction factor), wc in kcf, f'c in ksi
    # Note: 33000 * (0.150)^1.5 = 33000 * 0.05809 = 1917 (standard approximation)
    Ec = 33000.0 * 1.0 * (wc ** 1.5) * math.sqrt(fc)  # ksi

    # AASHTO LRFD Eq. 5.4.2.6-1: fr = 0.24 * sqrt(f'c) (ksi)
    fr = 0.24 * math.sqrt(fc)

    # ACI 318 / Hognestad: eps_0 = 2*f'c / Ec
    eps_0 = 2.0 * fc

def unconfined_concrete(tag: int, fc: float, eps_0: Optional[float] = None,
                        eps_cu: Optional[float] = None) -> int:
    """Define unconfined (cover) concrete using Concrete01.

    OpenSees Concrete01: zero tensile strength, linear-to-peak then linear
    degradation to crushing.

    Args:
        tag:    Material tag.
        fc:     Compressive strength f'c (ksi, positive). Internally negated
                for OpenSees convention.
        eps_0:  Strain at peak stress (positive). Default: computed from fc.
        eps_cu: Ultimate crushing strain (positive). Default: 0.003.

    Returns:
        Material tag.

    Reference:
        AASHTO LRFD 5.6.3.3.2 — strain limits for concrete.
    """
    props = concrete_defaults(fc)
    e0 = eps_0 if eps_0 is not None else props.eps_0
    ecu = eps_cu if eps_cu is not None else props.eps_cu

    # Concrete01: fpc, epsc0, fpcu, epsU
    # fpc = peak compressive stress (negative in OpenSees)
    # epsc0 = strain at p

def confined_concrete(tag: int, fc: float, fcc: float, ecc: float,
                      ecu: float) -> int:
    """Define confined concrete using Concrete01.

    Mander model: confined concrete has enhanced strength (fcc > fc) and
    enhanced ductility (ecu >> 0.003).

    Args:
        tag:  Material tag.
        fc:   Unconfined compressive strength (ksi, positive).
        fcc:  Confined compressive strength (ksi, positive).
        ecc:  Strain at confined peak stress (positive).
        ecu:  Ultimate confined crushing strain (positive).

    Returns:
        Material tag.

    Reference:
        Mander et al. (1988), AASHTO SGS Guide Spec 8.4.4.
    """
    # Residual strength: typically 0.2*fcc for confined concrete
    fcu_residual = 0.2 * fcc
    ops.uniaxialMaterial('Concrete01', tag, -fcc, -ecc, -fcu_residual, -ecu)
    return tag

def mander_confinement(fc: float, fy_transverse: float, rho_s: float,
                       config: str = "circular") -> Tuple[float, float, float]:
    """Compute confined concrete parameters per Mander et al. (1988).

    Args:
        fc:             Unconfined f'c (ksi, positive).
        fy_transverse:  Yield strength of transverse reinforcement (ksi).
        rho_s:          Volumetric ratio of transverse steel.
        config:         "circular" for spiral/hoop, "rectangular" for ties.

    Returns:
        Tuple of (fcc, ecc, ecu):
            fcc — confined compressive strength (ksi)
            ecc — strain at confined peak
            ecu — ultimate confined strain

    Reference:
        Mander, Priestley & Park (1988).
        AASHTO SGS Guide Spec Section 8.4.4.
    """
    # Effective confinement pressure
    # ke = confinement effectiveness coefficient
    if config == "circular":
        ke = 0.95  # spirals
    else:
        ke = 0.75  # rectangular ties (approximate

class SteelDefaults:
    """Common steel material property sets.

    All values in ksi.
    """
    fy: float       # Yield strength
    fu: float       # Ultimate strength
    Es: float       # Elastic modulus
    b: float        # Strain-hardening ratio (Esh/Es)
    R0: float = 20.0    # Steel02 transition parameter
    cR1: float = 0.925   # Steel02 transition parameter
    cR2: float = 0.15    # Steel02 transition parameter

def reinforcing_steel(tag: int, fy: float = 60.0, fu: float = 90.0,
                      Es: float = 29000.0, b: float = 0.01,
                      R0: float = 20.0, cR1: float = 0.925,
                      cR2: float = 0.15) -> int:
    """Define reinforcing steel using Steel02 (Giuffré-Menegotto-Pinto).

    Suitable for A615 Gr60 and A706 Gr60 reinforcement. Steel02 captures
    the Bauschinger effect under cyclic loading, critical for seismic analysis.

    Args:
        tag:  Material tag.
        fy:   Yield strength (ksi). Default: 60 ksi (Gr60).
        fu:   Ultimate strength (ksi). Default: 90 ksi.
        Es:   Elastic modulus (ksi). Default: 29000 ksi.
        b:    Strain-hardening ratio (b = Esh/Es). Default: 0.01.
        R0:   Initial value of curvature parameter. Default: 20.
        cR1:  Curvature degradation parameter. Default: 0.925.
        cR2:  Curvature degradation parameter. Default: 0.15.

    Returns:
        Material tag.

    Reference:
        AASHTO L

def structural_steel(tag: int, fy: float = 50.0, Es: float = 29000.0,
                     b: float = 0.01) -> int:
    """Define structural steel for W-shapes using Steel02.

    Default: ASTM A992 Grade 50 for rolled W-shapes.

    Args:
        tag: Material tag.
        fy:  Yield strength (ksi). Default: 50 ksi (A992).
        Es:  Elastic modulus (ksi). Default: 29000 ksi.
        b:   Strain-hardening ratio. Default: 0.01.

    Returns:
        Material tag.

    Reference:
        AASHTO LRFD 6.4.1: Structural Steel.
        ASTM A992/A992M — Standard Specification for Structural Steel Shapes.
    """
    ops.uniaxialMaterial('Steel02', tag, fy, Es, b, 20.0, 0.925, 0.15)
    return tag

def prestressing_strand(tag: int, fpu: float = 270.0,
                        Eps: float = 28500.0) -> int:
    """Define prestressing strand using Steel02.

    Calibrated for 7-wire, low-relaxation strand per ASTM A416.
    fpy = 0.9 * fpu for low-relaxation strand.

    Args:
        tag:  Material tag.
        fpu:  Ultimate tensile strength (ksi). Default: 270 ksi.
        Eps:  Elastic modulus (ksi). Default: 28500 ksi.

    Returns:
        Material tag.

    Reference:
        AASHTO LRFD 5.4.4: Prestressing Steel.
        ASTM A416/A416M: Standard Specification for Low-Relaxation,
        Seven-Wire Steel Strand for Prestressed Concrete.
    """
    fpy = 0.9 * fpu  # Yield strength for low-relaxation strand
    b = 0.008  # Low strain-hardening ratio for strand
    # R0=10 gives better fit for strand stress-strain curve
    ops.uniaxialMaterial('Steel02', tag, fpy, Eps, b, 10.0, 0.9, 0.15)
    return tag

class SoilLayer:
    """Soil layer definition for spring parameter computation.

    Attributes:
        depth_top:    Depth to top of layer (inches from mudline).
        depth_bot:    Depth to bottom of layer (inches).
        soil_type:    "sand" or "clay".
        gamma:        Effective unit weight (kip/in³). Typical sand: 0.0000694 (120 pcf).
        phi:          Friction angle (degrees). For sand.
        su:           Undrained shear strength (ksi). For clay.
        eps_50:       Strain at 50% of ultimate resistance. For clay.
        k_py:         Initial modulus of subgrade reaction (kip/in³). For sand.
    """
    depth_top: float
    depth_bot: float
    soil_type: str  # "sand" or "clay"
    gamma: float = 0.0000694  # ~120 pcf in kip/in³
    phi: float = 35.0
    su: float = 0.0
    eps_50: float = 0.01
    k_py: float = 0.0

def py_spring(tag: int, soil_type: int, pu: float, y50: float,
              cd: float = 0.0) -> int:
    """Define lateral soil spring using PySimple1.

    Models the p-y behavior for laterally loaded piles.

    Args:
        tag:        Material tag.
        soil_type:  1 = soft clay (Matlock 1970), 2 = sand (API).
        pu:         Ultimate capacity of the p-y spring (kip/in).
        y50:        Displacement at 50% of pu (inches).
        cd:         Drag resistance ratio (0 to 1). Default: 0.

    Returns:
        Material tag.

    Reference:
        API RP 2GEO Section 8.
        Boulanger, R.W. et al. (1999). "Seismic Soil-Pile-Structure
        Interaction Experiments and Analyses." ASCE JGGE.
    """
    ops.uniaxialMaterial('PySimple1', tag, soil_type, pu, y50, cd)
    return tag

def tz_spring(tag: int, soil_type: int, tult: float, z50: float) -> int:
    """Define skin friction spring using TzSimple1.

    Models the t-z behavior for axially loaded piles (shaft friction).

    Args:
        tag:        Material tag.
        soil_type:  1 = Reese & O'Neill t-z for driven piles,
                    2 = Mosher (1984) t-z for drilled shafts.
        tult:       Ultimate skin friction capacity (kip/in).
        z50:        Displacement at 50% of tult (inches).

    Returns:
        Material tag.

    Reference:
        API RP 2GEO Section 9.
        FHWA-NHI-10-016 Section 9.3: Axial Capacity.
    """
    ops.uniaxialMaterial('TzSimple1', tag, soil_type, tult, z50)
    return tag

def qz_spring(tag: int, qult: float, z50: float, soil_type: int = 2) -> int:
    """Define tip bearing spring using QzSimple1.

    Models the q-z behavior for pile tip resistance.

    Args:
        tag:        Material tag.
        qult:       Ultimate tip bearing capacity (kip).
        z50:        Displacement at 50% of qult (inches).
        soil_type:  1 = Reese & O'Neill backbone, 2 = Vijayvergiya (1977).

    Returns:
        Material tag.

    Reference:
        API RP 2GEO Section 9.
        FHWA-NHI-10-016: Drilled Shafts Chapter 13.
    """
    ops.uniaxialMaterial('QzSimple1', tag, soil_type, qult, z50)
    return tag

def api_py_curves(diameter: float, depth: float,
                  soil_layers: List[SoilLayer]) -> List[Dict[str, float]]:
    """Generate p-y curve parameters per API RP 2GEO / AASHTO.

    Computes pu and y50 at a given depth for the soil profile.

    Args:
        diameter:     Pile diameter (inches).
        depth:        Depth below mudline (inches).
        soil_layers:  List of SoilLayer definitions.

    Returns:
        List of dicts with keys: 'depth', 'pu', 'y50', 'soil_type_code'.

    Reference:
        API RP 2GEO (2014) Section 8.6-8.7.
        AASHTO LRFD Section 10.7.3.
    """
    results = []
    for layer in soil_layers:
        if depth < layer.depth_top or depth > layer.depth_bot:
            continue

        if layer.soil_type == "sand":
            pu, y50 = _api_py_sand(diameter, depth, layer)
            soil_code = 2
        else:  # clay
            pu, y50 = _api_py_clay(diameter, depth, layer)
            soil_code = 1

        results.append({
        

def _api_py_sand(diameter: float, depth: float,
                 layer: SoilLayer) -> Tuple[float, float]:
    """Compute p-y parameters for sand per API RP 2GEO Section 8.7.

    Args:
        diameter: Pile diameter (inches).
        depth:    Depth below mudline (inches).
        layer:    SoilLayer with sand properties.

    Returns:
        Tuple of (pu, y50) in kip/in and inches.
    """
    phi_rad = math.radians(layer.phi)
    gamma = layer.gamma  # kip/in³

    # Overburden pressure at depth
    sigma_v = gamma * depth  # ksi

    # API coefficients (simplified)
    # C1, C2, C3 depend on phi — use approximate relations
    c1 = max(0.0, (0.115 * 10 ** (0.0405 * layer.phi)))
    c2 = max(0.0, (0.571 * 10 ** (0.022 * layer.phi)))
    c3 = max(0.0, (0.646 * 10 ** (0.0555 * layer.phi)))

    # Shallow vs deep mechanism
    pu_shallow = (c1 * depth + c2 * diameter) * sigma_v
    pu_deep = c3 * sigma_v * diameter
    pu = min(pu_shallow, pu_deep) if depth > 0 else 0.001

    pu = m

def _api_py_clay(diameter: float, depth: float,
                 layer: SoilLayer) -> Tuple[float, float]:
    """Compute p-y parameters for soft clay per Matlock (1970).

    Args:
        diameter: Pile diameter (inches).
        depth:    Depth below mudline (inches).
        layer:    SoilLayer with clay properties (su, eps_50).

    Returns:
        Tuple of (pu, y50) in kip/in and inches.

    Reference:
        Matlock, H. (1970). "Correlations for Design of Laterally Loaded
        Piles in Soft Clay." OTC 1204.
    """
    su = layer.su  # ksi
    gamma = layer.gamma  # kip/in³
    J = 0.5  # Matlock's J factor (0.25-0.5, typical 0.5)

    # Ultimate soil resistance
    # Shallow: Np = 3 + gamma*z/su + J*z/D
    # Deep: Np = 9
    if su > 0 and diameter > 0:
        Np_shallow = 3.0 + (gamma * depth) / su + J * depth / diameter
    else:
        Np_shallow = 3.0
    Np = min(Np_shallow, 9.0)

    pu = Np * su * diameter  # kip/in (force per unit length)
    pu = max(pu, 0.001)

def _api_k_sand(phi: float) -> float:
    """Approximate initial modulus of subgrade reaction for sand (kip/in³).

    From API RP 2GEO Table 8.7-1 (interpolated).

    Args:
        phi: Friction angle in degrees.

    Returns:
        k in kip/in³.
    """
    # Approximate linear interpolation of API table
    # phi: 25→5, 30→10, 35→25, 40→50 (lb/in³ then /1000 for kip/in³)
    # Actually the API table gives k in lb/in³: 25→20, 30→42, 35→70, 40→140
    # Convert to kip/in³
    if phi <= 25:
        k_pci = 20.0
    elif phi <= 30:
        k_pci = 20.0 + (42.0 - 20.0) * (phi - 25.0) / 5.0
    elif phi <= 35:
        k_pci = 42.0 + (70.0 - 42.0) * (phi - 30.0) / 5.0
    elif phi <= 40:
        k_pci = 70.0 + (140.0 - 70.0) * (phi - 35.0) / 5.0
    else:
        k_pci = 140.0

    return k_pci / 1000.0  # Convert lb/in³ to kip/in³

def compute_tult(diameter: float, depth: float,
                 soil_layer: SoilLayer) -> float:
    """Compute ultimate unit skin friction (tult) for axially loaded piles.

    For clay: alpha method (FHWA-NHI-10-016 Section 9.3.2).
    For sand: beta method (FHWA-NHI-10-016 Section 9.3.3).

    Args:
        diameter:    Pile diameter (inches).
        depth:       Depth to midpoint of layer (inches).
        soil_layer:  SoilLayer definition.

    Returns:
        Ultimate unit skin friction tult (kip/in) — force per unit length
        of pile. Multiply by tributary length for total spring capacity.

    Reference:
        FHWA-NHI-10-016: Drilled Shafts, Chapter 9.
        AASHTO LRFD Section 10.8.3.5.
    """
    perimeter = math.pi * diameter

    if soil_layer.soil_type == "clay":
        # Alpha method: f_s = alpha * su
        su = soil_layer.su  # ksi
        # Alpha from O'Neill & Reese (1999):
        # su/pa <= 1.5: alpha = 0.55
        # su/pa > 1.5: alpha = 0.55 - 0.1*