Fix R degradation in SteelMPF when strain history starts at zero

[gaaraujo]

Fix R degradation in SteelMPF when strain history starts at zero

This PR fixes a bug in SteelMPF that caused incorrect degradation of the curvature parameter R when the strain history starts at $\varepsilon = 0$.

The issue led to inconsistent responses for equivalent loading histories depending on whether -prependZero was used in the Path time series.

Adapted from OpenSees #1734 / 97bccae.

Reproduction

Current behavior	Behavior with this PR

Test script

import matplotlib.pyplot as plt
import numpy as np
import xara

E = 29000.0
FY = 50.0
B_P = 0.02
B_N = 0.02
R0 = 20.0
C_R1 = 0.925
C_R2 = 0.15
EPS_Y = FY / E

EPS_PEAK = 2.0 * EPS_Y
N = 10
D = np.linspace(0.0, EPS_PEAK, N, endpoint=False)
DISP_HISTORY = np.r_[D, EPS_PEAK - D, -D, -(EPS_PEAK - D)]


def _scalar(response):
    return float(response[0] if hasattr(response, "__getitem__") else response)


def run_simulation(disp_history, prepend_zero=False):
    model = xara.Model(ndm=1, ndf=1)
    model.node(1, (0.0,))
    model.fix(1, (1,))
    model.node(2, (0.0,))
    model.uniaxialMaterial(
        "SteelMPF", 1,
        FY, FY, E, B_P, B_N, R0, C_R1, C_R2,
    )
    model.element("zeroLength", 1, (1, 2), mat=1, dir=1)

    dt = 1.0
    extras = {"useLast": True}
    if prepend_zero:
        extras["prependZero"] = True
    model.timeSeries("Path", 1, dt=dt, values=list(disp_history), **extras)
    model.pattern("Plain", 1, 1)
    model.sp(2, 1, 1.0, pattern=1)
    model.integrator("LoadControl", dt)
    model.constraints("Transformation")
    model.numberer("Plain")
    model.system("UmfPack")
    model.analysis("Static", "-noWarnings")

    n = len(disp_history)
    stress = np.zeros(n)
    strain = np.zeros(n)
    time = np.zeros(n)
    for i in range(n):
        if model.analyze(1) != 0:
            raise RuntimeError("analyze failed at step %d/%d" % (i + 1, n))
        stress[i] = _scalar(model.eleResponse(1, "material", 1, "stress"))
        strain[i] = _scalar(model.eleResponse(1, "material", 1, "strain"))
        time[i] = model.getTime()
    return stress, strain, time


def main():
    runs = (
        (run_simulation(DISP_HISTORY, False), "prepend_zero=False", "-x"),
        (run_simulation(DISP_HISTORY, True), "prepend_zero=True", "--o"),
    )

    fig, (ax_t, ax_ss) = plt.subplots(2, 1, figsize=(5.5, 6.0))
    for (s, e, t), lab, sty in runs:
        en = e / EPS_Y
        ax_t.plot(t, en, sty, label=lab, lw=1.2, markersize=4)
        ax_ss.plot(en, s / FY, sty, label=lab, lw=1.2, markersize=4)

    ax_t.set_ylabel(r"$\varepsilon / \varepsilon_y$")
    ax_t.set_xlabel("pseudo-time")
    ax_t.axhline(0.0, color="0.8", lw=0.8)
    ax_t.grid(True, alpha=0.35)

    ax_ss.set_xlabel(r"$\varepsilon / \varepsilon_y$")
    ax_ss.set_ylabel(r"$\sigma / f_y$")
    ax_ss.axhline(0.0, color="0.8", lw=0.8)
    ax_ss.axvline(0.0, color="0.8", lw=0.8)
    ax_ss.grid(True, alpha=0.35)
    ax_ss.legend(loc="best", fontsize=8)

    fig.tight_layout()
    fig.savefig("demo_prepend_zero_bug.png", dpi=150)
    plt.close(fig)
    print("Wrote demo_prepend_zero_bug.png")


if __name__ == "__main__":
    main()