"""Integration: verify zeroLengthSection runs and traces moment-curvature. Mirrors OpenSees's Example 2 (Moment-Curvature of a rectangular RC section) but with a simplified elastic material so we can check the slope against a closed-form value. Full Concrete01/Steel01 fiber behaviour is exercised by the Phase 9 pushover tests. """ from __future__ import annotations import pytest pytest.importorskip("openseespy") from otko.core import ( ElasticUniaxial, FiberSection, LinearTimeSeries, NodalLoad, Node, PlainLoadPattern, Project, RectangularPatch, StaticCase, ZeroLengthSectionElement, ) from otko.services.opensees_runner import OpenSeesRunner pytestmark = pytest.mark.slow def _moment_curvature_project(moment: float) -> Project: """Two coincident nodes + a rectangular fibre section + one moment step. Node 1 clamped; node 2 free in Ux and Rz. Applied moment = ``moment`` at node 2's DOF 3 (Rz). With a linear-elastic fibre material the curvature should be ``moment / (E·I)``. """ E = 30000.0 # Elastic modulus b, h = 10.0, 20.0 # width × depth (in) return Project( ndm=2, ndf=3, nodes=[ Node(id=1, coords=(0, 0, 0), restraint=(True, True, False, False, False, True)), Node(id=2, coords=(0, 0, 0), restraint=(False, True, False, False, False, False)), ], materials=[ElasticUniaxial(id=1, name="Elastic", E=E)], sections=[ FiberSection( id=1, name="Rect", patches=[ RectangularPatch( material_id=1, n_fib_y=20, n_fib_z=1, y_i=-h / 2, z_i=-b / 2, y_j=h / 2, z_j=b / 2, ) ], ) ], elements=[ZeroLengthSectionElement(id=1, nodes=(1, 2), section_id=1)], time_series=[LinearTimeSeries(id=1, name="R")], load_patterns=[ PlainLoadPattern( id=1, time_series_id=1, # NodalLoad.forces = (Fx, Fy, Fz, Mx, My, Mz). Moment around # z (= curvature driver in 2D) goes into index 5, not 2. nodal_loads=[NodalLoad(node_id=2, forces=(0, 0, 0, 0, 0, moment))], ) ], analyses=[StaticCase(id=1, name="MK", pattern_ids=[1], n_steps=1)], ) def test_zero_length_section_elastic_curvature_matches_closed_form() -> None: """For a linear fibre section, curvature = M / (E·I).""" # Moment M → curvature κ = M / (E·I). Rectangle: I = b·h³/12. M = 500.0 E = 30000.0 b, h = 10.0, 20.0 I = b * h**3 / 12.0 expected_kappa = M / (E * I) proj = _moment_curvature_project(moment=M) result = OpenSeesRunner(proj).run(proj.analyses[0]) # Rz at node 2 IS the curvature for a zero-length section. ux, uy, rz = result.node_disp[2][-1] assert rz == pytest.approx( expected_kappa, rel=5e-3 ), f"κ = {rz:.6e}, expected {expected_kappa:.6e}" def test_pushover_drives_rotation_for_moment_curvature() -> None: """Full moment-curvature analysis via PushoverCase on DOF 3 (Rz). Mirrors the OpenSees Moment-Curvature example's driver: a DisplacementControl pushover on node 2's rotational DOF produces a moment-curvature curve. For a linear-elastic fibre section the base "shear" is actually the reactive moment, and the curve is a straight line through the origin with slope E·I. """ from otko.core import PushoverCase E = 30000.0 b, h = 10.0, 20.0 I = b * h**3 / 12.0 target_kappa = 1e-5 steps = 20 # PushoverCase with DisplacementControl scales the load pattern — # needs a *non-zero* reference moment at the control DOF. proj = _moment_curvature_project(moment=1.0) proj.analyses = [ PushoverCase( id=1, name="MK-push", pattern_ids=[1], control_node=2, control_dof=3, # DOF 3 = Rz target_disp=target_kappa, # "displacement" == curvature here step_size=target_kappa / steps, base_nodes=[1], ) ] result = OpenSeesRunner(proj).run(proj.analyses[0]) # Every (κ, M) point must satisfy M = E·I·κ (1 % tolerance allows # for the ~20-fibre discretisation of the rectangular section). for kappa, moment in zip(result.control_disp, result.base_shear, strict=False): if abs(kappa) < 1e-12: continue expected_M = E * I * kappa assert moment == pytest.approx( expected_M, rel=1e-2 ), f"at κ={kappa:.3e}: M={moment:.3e}, expected {expected_M:.3e}" # Terminal curvature must reach the target. assert result.control_disp[-1] == pytest.approx(target_kappa, rel=1e-3) def test_moment_curvature_with_constant_axial_preload() -> None: """OpenSees MK Example 2 recipe: Concrete01 + Steel01 fibre section, constant axial compression preloaded, then DisplacementControl ramps curvature. Verifies the runner's two-stage preload + pushover plumbing (the key fix that makes convergence possible on nonlinear RC sections). """ from otko.core import ( Concrete01, ConstantTimeSeries, FiberSection, LinearTimeSeries, NodalLoad, PlainLoadPattern, PushoverCase, RectangularPatch, Steel01, StraightLayer, ) colWidth = 15.0 colDepth = 24.0 cover = 1.5 As = 0.60 y1 = colDepth / 2 z1 = colWidth / 2 proj = Project( ndm=2, ndf=3, nodes=[ Node(id=1, coords=(0, 0, 0), restraint=(True, True, False, False, False, True)), Node(id=2, coords=(0, 0, 0), restraint=(False, True, False, False, False, False)), ], materials=[ Concrete01(id=1, name="Core", fpc=-6.0, epsc0=-0.004, fpcu=-5.0, epsU=-0.014), Concrete01(id=2, name="Cover", fpc=-5.0, epsc0=-0.002, fpcu=0.0, epsU=-0.006), Steel01(id=3, name="Steel", Fy=60.0, E0=30000.0, b=0.01), ], sections=[ FiberSection( id=1, name="RC", patches=[ # Core (confined) RectangularPatch( material_id=1, n_fib_y=10, n_fib_z=1, y_i=cover - y1, z_i=cover - z1, y_j=y1 - cover, z_j=z1 - cover, ), # Top cover RectangularPatch( material_id=2, n_fib_y=10, n_fib_z=1, y_i=-y1, z_i=z1 - cover, y_j=y1, z_j=z1, ), # Bottom cover RectangularPatch( material_id=2, n_fib_y=10, n_fib_z=1, y_i=-y1, z_i=-z1, y_j=y1, z_j=cover - z1, ), # Left cover RectangularPatch( material_id=2, n_fib_y=2, n_fib_z=1, y_i=-y1, z_i=cover - z1, y_j=cover - y1, z_j=z1 - cover, ), # Right cover RectangularPatch( material_id=2, n_fib_y=2, n_fib_z=1, y_i=y1 - cover, z_i=cover - z1, y_j=y1, z_j=z1 - cover, ), ], layers=[ StraightLayer( material_id=3, n_bars=3, bar_area=As, y_start=y1 - cover, z_start=z1 - cover, y_end=y1 - cover, z_end=cover - z1, ), StraightLayer( material_id=3, n_bars=2, bar_area=As, y_start=0.0, z_start=z1 - cover, y_end=0.0, z_end=cover - z1, ), StraightLayer( material_id=3, n_bars=3, bar_area=As, y_start=cover - y1, z_start=z1 - cover, y_end=cover - y1, z_end=cover - z1, ), ], ) ], elements=[ZeroLengthSectionElement(id=1, nodes=(1, 2), section_id=1)], time_series=[ ConstantTimeSeries(id=1, name="AxialP"), LinearTimeSeries(id=2, name="RefMoment"), ], load_patterns=[ PlainLoadPattern( id=1, name="AxialP", time_series_id=1, nodal_loads=[NodalLoad(node_id=2, forces=(-180.0, 0, 0, 0, 0, 0))], ), PlainLoadPattern( id=2, name="RefMoment", time_series_id=2, nodal_loads=[NodalLoad(node_id=2, forces=(0, 0, 0, 0, 0, 1.0))], ), ], analyses=[], ) # Yield curvature estimate from the Tcl example. d = colDepth - cover Ky = 60.0 / 30000.0 / (0.7 * d) target = Ky * 15 # μ = 15 proj.analyses = [ PushoverCase( id=1, name="MK", pattern_ids=[1, 2], control_node=2, control_dof=3, target_disp=target, step_size=target / 100, base_nodes=[1], test="NormUnbalance", tolerance=1e-9, max_iter=25, ) ] result = OpenSeesRunner(proj).run(proj.analyses[0]) # Analysis must actually converge past yield (not collapse at # step 1 like it did before the two-stage preload fix). assert len(result.control_disp) > 50, ( f"Converged for only {len(result.control_disp)} of 100 steps — " "preload stage broken?" ) # Curvature reached or passed yield. kappa_max = float(max(abs(k) for k in result.control_disp)) assert kappa_max > Ky, f"κ_max={kappa_max:.3e} < Ky={Ky:.3e}" # Nonlinear → curve has a distinct softening: slope late in the # run should be smaller than slope near the origin. d_early = (result.base_shear[5] - result.base_shear[1]) / ( result.control_disp[5] - result.control_disp[1] ) last = len(result.control_disp) - 1 mid = last // 2 d_late = (result.base_shear[last] - result.base_shear[mid]) / ( result.control_disp[last] - result.control_disp[mid] ) assert abs(d_late) < abs(d_early), ( f"Late slope {d_late:.2e} not smaller than early {d_early:.2e} " "— section response looks linear, preload probably didn't apply." )