carapace_core/model/materials/
hysteretic.rs

1use super::Material;
2
3/// `HystereticMaterial`'s `CloadIndicator`/`TloadIndicator` — which
4/// direction the material last loaded in (`None` only before the first
5/// nonzero strain increment ever).
6#[derive(Debug, Clone, Copy, PartialEq, Eq)]
7pub enum HystereticLoad {
8    None,
9    Positive,
10    Negative,
11}
12
13/// The full field set behind `Material::Hysteretic` (see that variant's
14/// doc comment for why it's boxed rather than inlined). `Copy` despite the
15/// field count — every field is itself `f64`/`HystereticLoad`, so this is
16/// cheap to copy out of its `Box` for the duration of one `evaluate` call,
17/// which is simpler than borrow-juggling a destructure through `&Box<_>`.
18///
19/// The backbone (`mom1p..mom3n`/`rot1p..rot3n`) and pinching/damage
20/// parameters (`pinch_x`, `pinch_y`, `damfc1`, `damfc2`, `beta`) are fixed
21/// at construction, along with the envelope slopes/limit (`e1p..e3n`,
22/// `eup`, `eun`, `energy_a`) `HystereticMaterial::setEnvelope` derives from
23/// them once — recomputing those every `evaluate` call would be wasted
24/// work, since they never change. History: `rot_max`/`rot_min` are the
25/// peak strain reached on each side (`CrotMax`/`CrotMin`); `rot_pu`/
26/// `rot_nu` are the pinching reference points captured at the last
27/// reversal into negative/positive loading respectively (`CrotPu`/
28/// `CrotNu`); `energy_d` is cumulative dissipated energy, used only to
29/// compute the damage factor on a reversal; `load` mirrors
30/// `CloadIndicator`.
31#[derive(Debug, Clone, Copy)]
32pub struct HystereticFields {
33    pub mom1p: f64,
34    pub rot1p: f64,
35    pub mom2p: f64,
36    pub rot2p: f64,
37    pub mom3p: f64,
38    pub rot3p: f64,
39    pub mom1n: f64,
40    pub rot1n: f64,
41    pub mom2n: f64,
42    pub rot2n: f64,
43    pub mom3n: f64,
44    pub rot3n: f64,
45    pub pinch_x: f64,
46    pub pinch_y: f64,
47    pub damfc1: f64,
48    pub damfc2: f64,
49    pub beta: f64,
50    pub e1p: f64,
51    pub e2p: f64,
52    pub e3p: f64,
53    pub e1n: f64,
54    pub e2n: f64,
55    pub e3n: f64,
56    pub eup: f64,
57    pub eun: f64,
58    pub energy_a: f64,
59    pub rot_max: f64,
60    pub rot_min: f64,
61    pub rot_pu: f64,
62    pub rot_nu: f64,
63    pub energy_d: f64,
64    pub load: HystereticLoad,
65    pub strain: f64,
66    pub stress: f64,
67    /// Committed tangent (OpenSees' `Ttangent` as left by the last trial),
68    /// returned unchanged for a trial at the committed strain.
69    pub tangent: f64,
70}
71
72/// `HystereticMaterial`'s fixed (non-history) backbone parameters, bundled
73/// so `hysteretic_positive_increment`/`hysteretic_negative_increment` don't
74/// need a 20-parameter signature each. Mirrors `posEnvlpStress`/
75/// `negEnvlpStress`/`posEnvlpTangent`/`negEnvlpTangent`/`posEnvlpRotlim`/
76/// `negEnvlpRotlim` from `HystereticMaterial.cpp` exactly.
77struct HystereticEnvelope {
78    mom1p: f64,
79    rot1p: f64,
80    mom2p: f64,
81    rot2p: f64,
82    mom3p: f64,
83    rot3p: f64,
84    e1p: f64,
85    e2p: f64,
86    e3p: f64,
87    mom1n: f64,
88    rot1n: f64,
89    mom2n: f64,
90    rot2n: f64,
91    mom3n: f64,
92    rot3n: f64,
93    e1n: f64,
94    e2n: f64,
95    e3n: f64,
96}
97
98/// `HystereticMaterial.cpp`'s `POS_INF_STRAIN`/`NEG_INF_STRAIN` sentinels
99/// (from `UniaxialMaterial.h`) — a strain limit this large means "no
100/// limit," never an actual physical bound.
101const HYSTERETIC_POS_INF_STRAIN: f64 = 1.0e16;
102const HYSTERETIC_NEG_INF_STRAIN: f64 = -1.0e16;
103
104impl HystereticEnvelope {
105    fn pos_stress(&self, strain: f64) -> f64 {
106        if strain <= 0.0 {
107            0.0
108        } else if strain <= self.rot1p {
109            self.e1p * strain
110        } else if strain <= self.rot2p {
111            self.mom1p + self.e2p * (strain - self.rot1p)
112        } else if strain <= self.rot3p || self.e3p > 0.0 {
113            self.mom2p + self.e3p * (strain - self.rot2p)
114        } else {
115            self.mom3p
116        }
117    }
118
119    fn neg_stress(&self, strain: f64) -> f64 {
120        if strain >= 0.0 {
121            0.0
122        } else if strain >= self.rot1n {
123            self.e1n * strain
124        } else if strain >= self.rot2n {
125            self.mom1n + self.e2n * (strain - self.rot1n)
126        } else if strain >= self.rot3n || self.e3n > 0.0 {
127            self.mom2n + self.e3n * (strain - self.rot2n)
128        } else {
129            self.mom3n
130        }
131    }
132
133    fn pos_tangent(&self, strain: f64) -> f64 {
134        if strain < 0.0 {
135            self.e1p * 1.0e-9
136        } else if strain <= self.rot1p {
137            self.e1p
138        } else if strain <= self.rot2p {
139            self.e2p
140        } else if strain <= self.rot3p || self.e3p > 0.0 {
141            self.e3p
142        } else {
143            self.e1p * 1.0e-9
144        }
145    }
146
147    fn neg_tangent(&self, strain: f64) -> f64 {
148        if strain > 0.0 {
149            self.e1n * 1.0e-9
150        } else if strain >= self.rot1n {
151            self.e1n
152        } else if strain >= self.rot2n {
153            self.e2n
154        } else if strain >= self.rot3n || self.e3n > 0.0 {
155            self.e3n
156        } else {
157            self.e1n * 1.0e-9
158        }
159    }
160
161    /// The strain at which the positive backbone, extended past `strain`,
162    /// would cross zero stress on a softening segment — `POS_INF_STRAIN`
163    /// if the backbone never softens (or hasn't yet, from `strain`).
164    ///
165    /// The two branches below returning the same `HYSTERETIC_POS_INF_STRAIN`
166    /// literal are checking genuinely different conditions from
167    /// `HystereticMaterial.cpp`'s `posEnvlpRotlim` (no strain limit was
168    /// ever found vs. the one found doesn't correspond to a real softening
169    /// crossing) — collapsing them would lose that distinction.
170    #[allow(clippy::if_same_then_else)]
171    fn pos_rotlim(&self, strain: f64) -> f64 {
172        let mut strain_limit = HYSTERETIC_POS_INF_STRAIN;
173
174        if strain <= self.rot1p {
175            return HYSTERETIC_POS_INF_STRAIN;
176        }
177        if strain > self.rot1p && strain <= self.rot2p && self.e2p < 0.0 {
178            strain_limit = self.rot1p - self.mom1p / self.e2p;
179        }
180        if strain > self.rot2p && self.e3p < 0.0 {
181            strain_limit = self.rot2p - self.mom2p / self.e3p;
182        }
183
184        if strain_limit == HYSTERETIC_POS_INF_STRAIN {
185            HYSTERETIC_POS_INF_STRAIN
186        } else if self.pos_stress(strain_limit) > 0.0 {
187            HYSTERETIC_POS_INF_STRAIN
188        } else {
189            strain_limit
190        }
191    }
192
193    /// See `pos_rotlim` above for why the same-looking branches aren't
194    /// collapsed.
195    #[allow(clippy::if_same_then_else)]
196    fn neg_rotlim(&self, strain: f64) -> f64 {
197        let mut strain_limit = HYSTERETIC_NEG_INF_STRAIN;
198
199        if strain >= self.rot1n {
200            return HYSTERETIC_NEG_INF_STRAIN;
201        }
202        if strain < self.rot1n && strain >= self.rot2n && self.e2n < 0.0 {
203            strain_limit = self.rot1n - self.mom1n / self.e2n;
204        }
205        if strain < self.rot2n && self.e3n < 0.0 {
206            strain_limit = self.rot2n - self.mom2n / self.e3n;
207        }
208
209        if strain_limit == HYSTERETIC_NEG_INF_STRAIN {
210            HYSTERETIC_NEG_INF_STRAIN
211        } else if self.neg_stress(strain_limit) < 0.0 {
212            HYSTERETIC_NEG_INF_STRAIN
213        } else {
214            strain_limit
215        }
216    }
217}
218
219impl Material {
220    /// The 3-point-per-side backbone constructor (OpenSees also offers a
221    /// 2-point constructor that derives the midpoint as the average of the
222    /// two given points — not ported, since a caller can just pass that
223    /// average as `mom2p`/`rot2p`/`mom2n`/`rot2n` directly).
224    ///
225    /// Debug builds assert the backbone is monotonic and one-to-one
226    /// (`0 < rot1p < rot2p < rot3p`, `rot1n < rot2n < rot3n < 0`) —
227    /// OpenSees hard-exits the process on this same precondition failure,
228    /// which isn't a library-appropriate response here.
229    #[allow(clippy::too_many_arguments)]
230    pub fn hysteretic(
231        mom1p: f64,
232        rot1p: f64,
233        mom2p: f64,
234        rot2p: f64,
235        mom3p: f64,
236        rot3p: f64,
237        mom1n: f64,
238        rot1n: f64,
239        mom2n: f64,
240        rot2n: f64,
241        mom3n: f64,
242        rot3n: f64,
243        pinch_x: f64,
244        pinch_y: f64,
245        damfc1: f64,
246        damfc2: f64,
247        beta: f64,
248    ) -> Self {
249        debug_assert!(
250            rot1p > 0.0 && rot2p > rot1p && rot3p > rot2p,
251            "positive backbone must be increasing"
252        );
253        debug_assert!(
254            rot1n < 0.0 && rot2n < rot1n && rot3n < rot2n,
255            "negative backbone must be decreasing"
256        );
257
258        let energy_a = 0.5
259            * (rot1p * mom1p
260                + (rot2p - rot1p) * (mom2p + mom1p)
261                + (rot3p - rot2p) * (mom3p + mom2p)
262                + rot1n * mom1n
263                + (rot2n - rot1n) * (mom2n + mom1n)
264                + (rot3n - rot2n) * (mom3n + mom2n));
265
266        let e1p = mom1p / rot1p;
267        let e2p = (mom2p - mom1p) / (rot2p - rot1p);
268        let e3p = (mom3p - mom2p) / (rot3p - rot2p);
269        let e1n = mom1n / rot1n;
270        let e2n = (mom2n - mom1n) / (rot2n - rot1n);
271        let e3n = (mom3n - mom2n) / (rot3n - rot2n);
272        let eup = e1p.max(e2p).max(e3p);
273        let eun = e1n.max(e2n).max(e3n);
274
275        Material::Hysteretic(Box::new(HystereticFields {
276            mom1p,
277            rot1p,
278            mom2p,
279            rot2p,
280            mom3p,
281            rot3p,
282            mom1n,
283            rot1n,
284            mom2n,
285            rot2n,
286            mom3n,
287            rot3n,
288            pinch_x,
289            pinch_y,
290            damfc1,
291            damfc2,
292            beta,
293            e1p,
294            e2p,
295            e3p,
296            e1n,
297            e2n,
298            e3n,
299            eup,
300            eun,
301            energy_a,
302            rot_max: 0.0,
303            rot_min: 0.0,
304            rot_pu: 0.0,
305            rot_nu: 0.0,
306            energy_d: 0.0,
307            load: HystereticLoad::None,
308            strain: 0.0,
309            stress: 0.0,
310            tangent: e1p,
311        }))
312    }
313}
314
315/// Ported from `HystereticMaterial::setTrialStrain`
316/// (`HystereticMaterial.cpp`). OpenSees' `TloadIndicator == 0 && strain ==
317/// 0.0` guard isn't ported: it returns the untouched initial
318/// `(Tstress, Ttangent) = (0.0, E1p)`, which the `|dStrain| < DBL_EPSILON`
319/// return below reproduces from the committed state (checked by
320/// `hysteretic_reproduces_initial_tangent_at_zero_strain`).
321pub(super) fn evaluate(m: &Material, strain: f64) -> (f64, f64, Material) {
322    let Material::Hysteretic(fields) = m else {
323        unreachable!()
324    };
325    let HystereticFields {
326        mom1p,
327        rot1p,
328        mom2p,
329        rot2p,
330        mom3p,
331        rot3p,
332        mom1n,
333        rot1n,
334        mom2n,
335        rot2n,
336        mom3n,
337        rot3n,
338        pinch_x,
339        pinch_y,
340        damfc1,
341        damfc2,
342        beta,
343        e1p,
344        e2p,
345        e3p,
346        e1n,
347        e2n,
348        e3n,
349        eup,
350        eun,
351        energy_a,
352        mut rot_max,
353        mut rot_min,
354        mut rot_pu,
355        mut rot_nu,
356        energy_d: energy_d_committed,
357        mut load,
358        strain: cstrain,
359        stress: cstress,
360        tangent: ctangent,
361    } = **fields;
362
363    let envelope = HystereticEnvelope {
364        mom1p,
365        rot1p,
366        mom2p,
367        rot2p,
368        mom3p,
369        rot3p,
370        e1p,
371        e2p,
372        e3p,
373        mom1n,
374        rot1n,
375        mom2n,
376        rot2n,
377        mom3n,
378        rot3n,
379        e1n,
380        e2n,
381        e3n,
382    };
383
384    let dstrain = strain - cstrain;
385
386    // `HystereticMaterial::setTrialStrain` returns early below
387    // `DBL_EPSILON`, leaving the previous `(Tstress, Ttangent)`.
388    if dstrain.abs() < f64::EPSILON {
389        return (cstress, ctangent, m.clone());
390    }
391
392    if load == HystereticLoad::None {
393        load = if dstrain < 0.0 {
394            HystereticLoad::Negative
395        } else {
396            HystereticLoad::Positive
397        };
398    }
399
400    let (stress, tangent) = if strain >= rot_max {
401        rot_max = strain;
402        load = HystereticLoad::Positive;
403        (envelope.pos_stress(strain), envelope.pos_tangent(strain))
404    } else if strain <= rot_min {
405        rot_min = strain;
406        load = HystereticLoad::Negative;
407        (envelope.neg_stress(strain), envelope.neg_tangent(strain))
408    } else if dstrain < 0.0 {
409        let (stress, tangent, new_rot_min, new_rot_pu, new_load) = hysteretic_negative_increment(
410            dstrain,
411            strain,
412            cstrain,
413            cstress,
414            rot_max,
415            rot_min,
416            rot_pu,
417            load,
418            rot1p,
419            rot1n,
420            beta,
421            eup,
422            eun,
423            pinch_x,
424            pinch_y,
425            damfc1,
426            damfc2,
427            energy_a,
428            energy_d_committed,
429            &envelope,
430        );
431        rot_min = new_rot_min;
432        rot_pu = new_rot_pu;
433        load = new_load;
434        (stress, tangent)
435    } else {
436        // `dstrain > 0.0` (an increment under `DBL_EPSILON` returned above).
437        let (stress, tangent, new_rot_max, new_rot_nu, new_load) = hysteretic_positive_increment(
438            dstrain,
439            strain,
440            cstrain,
441            cstress,
442            rot_max,
443            rot_min,
444            rot_nu,
445            load,
446            rot1p,
447            rot1n,
448            beta,
449            eup,
450            eun,
451            pinch_x,
452            pinch_y,
453            damfc1,
454            damfc2,
455            energy_a,
456            energy_d_committed,
457            &envelope,
458        );
459        rot_max = new_rot_max;
460        rot_nu = new_rot_nu;
461        load = new_load;
462        (stress, tangent)
463    };
464
465    let energy_d = energy_d_committed + 0.5 * (cstress + stress) * dstrain;
466
467    (
468        stress,
469        tangent,
470        Material::Hysteretic(Box::new(HystereticFields {
471            mom1p,
472            rot1p,
473            mom2p,
474            rot2p,
475            mom3p,
476            rot3p,
477            mom1n,
478            rot1n,
479            mom2n,
480            rot2n,
481            mom3n,
482            rot3n,
483            pinch_x,
484            pinch_y,
485            damfc1,
486            damfc2,
487            beta,
488            e1p,
489            e2p,
490            e3p,
491            e1n,
492            e2n,
493            e3n,
494            eup,
495            eun,
496            energy_a,
497            rot_max,
498            rot_min,
499            rot_pu,
500            rot_nu,
501            energy_d,
502            load,
503            strain,
504            stress,
505            tangent,
506        })),
507    )
508}
509
510/// Ported from `HystereticMaterial::positiveIncrement`
511/// (`HystereticMaterial.cpp`). Called only for the "interior" branch —
512/// `rot_min < strain < rot_max` hasn't been hit yet by the caller's
513/// envelope-extreme checks. `rot_max`/`rot_min`/`rot_nu`/`load` are the
514/// values *entering* this call (`rot_min` and `rot_max` are read but only
515/// `rot_max` is ever updated here); returns
516/// `(stress, tangent, new_rot_max, new_rot_nu, new_load)`.
517#[allow(clippy::too_many_arguments)]
518fn hysteretic_positive_increment(
519    dstrain: f64,
520    strain: f64,
521    cstrain: f64,
522    cstress: f64,
523    rot_max: f64,
524    rot_min: f64,
525    rot_nu: f64,
526    mut load: HystereticLoad,
527    rot1p: f64,
528    rot1n: f64,
529    beta: f64,
530    eup: f64,
531    eun: f64,
532    pinch_x: f64,
533    pinch_y: f64,
534    damfc1: f64,
535    damfc2: f64,
536    energy_a: f64,
537    energy_d: f64,
538    envelope: &HystereticEnvelope,
539) -> (f64, f64, f64, f64, HystereticLoad) {
540    let kn_raw = (rot_min / rot1n).powf(beta);
541    let kn = if kn_raw < 1.0 { 1.0 } else { 1.0 / kn_raw };
542    let kp_raw = (rot_max / rot1p).powf(beta);
543    let kp = if kp_raw < 1.0 { 1.0 } else { 1.0 / kp_raw };
544
545    let mut rot_max = rot_max;
546    let mut rot_nu = rot_nu;
547
548    // Reversal from negative to positive loading: capture the pinching
549    // reference point `rot_nu` and apply ductility-based damage to
550    // `rot_max`, but only if the material had actually crossed into
551    // tension-side stress (`cstress <= 0.0`) before reversing.
552    if load == HystereticLoad::Negative && cstress <= 0.0 {
553        rot_nu = cstrain - cstress / (eun * kn);
554        let energy = energy_d - 0.5 * cstress / (eun * kn) * cstress;
555        let mut damfc = 0.0;
556        if rot_min < rot1n {
557            damfc = damfc2 * energy / energy_a;
558            damfc += damfc1 * (rot_min - rot1n) / rot1n;
559        }
560        rot_max *= 1.0 + damfc;
561    }
562    load = HystereticLoad::Positive;
563
564    if rot_max > HYSTERETIC_POS_INF_STRAIN {
565        rot_max = HYSTERETIC_POS_INF_STRAIN;
566    }
567    rot_max = if rot_max > rot1p { rot_max } else { rot1p };
568
569    let maxmom = envelope.pos_stress(rot_max);
570    let rotlim = envelope.neg_rotlim(rot_min);
571    let rotrel = if rotlim > rot_nu { rotlim } else { rot_nu };
572
573    let rotmp2 = rot_max - (1.0 - pinch_y) * maxmom / (eup * kp);
574    let rotch = rotrel + (rotmp2 - rotrel) * pinch_x;
575
576    let (stress, tangent) = if strain < rot_nu {
577        let t = eun * kn;
578        let s = cstress + t * dstrain;
579        if s >= 0.0 {
580            (0.0, eun * 1.0e-9)
581        } else {
582            (s, t)
583        }
584    } else if strain < rotch {
585        if strain <= rotrel {
586            (0.0, eup * 1.0e-9)
587        } else {
588            let t = maxmom * pinch_y / (rotch - rotrel);
589            let tmpmo1 = cstress + eup * kp * dstrain;
590            let tmpmo2 = (strain - rotrel) * t;
591            if tmpmo1 < tmpmo2 {
592                (tmpmo1, eup * kp)
593            } else {
594                (tmpmo2, t)
595            }
596        }
597    } else {
598        let t = (1.0 - pinch_y) * maxmom / (rot_max - rotch);
599        let tmpmo1 = cstress + eup * kp * dstrain;
600        let tmpmo2 = pinch_y * maxmom + (strain - rotch) * t;
601        if tmpmo1 < tmpmo2 {
602            (tmpmo1, eup * kp)
603        } else {
604            (tmpmo2, t)
605        }
606    };
607
608    (stress, tangent, rot_max, rot_nu, load)
609}
610
611/// Ported from `HystereticMaterial::negativeIncrement`
612/// (`HystereticMaterial.cpp`) — the mirror image of
613/// `hysteretic_positive_increment` above; see its doc comment.
614#[allow(clippy::too_many_arguments)]
615fn hysteretic_negative_increment(
616    dstrain: f64,
617    strain: f64,
618    cstrain: f64,
619    cstress: f64,
620    rot_max: f64,
621    rot_min: f64,
622    rot_pu: f64,
623    mut load: HystereticLoad,
624    rot1p: f64,
625    rot1n: f64,
626    beta: f64,
627    eup: f64,
628    eun: f64,
629    pinch_x: f64,
630    pinch_y: f64,
631    damfc1: f64,
632    damfc2: f64,
633    energy_a: f64,
634    energy_d: f64,
635    envelope: &HystereticEnvelope,
636) -> (f64, f64, f64, f64, HystereticLoad) {
637    let kn_raw = (rot_min / rot1n).powf(beta);
638    let kn = if kn_raw < 1.0 { 1.0 } else { 1.0 / kn_raw };
639    let kp_raw = (rot_max / rot1p).powf(beta);
640    let kp = if kp_raw < 1.0 { 1.0 } else { 1.0 / kp_raw };
641
642    let mut rot_min = rot_min;
643    let mut rot_pu = rot_pu;
644
645    if load == HystereticLoad::Positive && cstress >= 0.0 {
646        rot_pu = cstrain - cstress / (eup * kp);
647        let energy = energy_d - 0.5 * cstress / (eup * kp) * cstress;
648        let mut damfc = 0.0;
649        if rot_max > rot1p {
650            damfc = damfc2 * energy / energy_a;
651            damfc += damfc1 * (rot_max - rot1p) / rot1p;
652        }
653        rot_min *= 1.0 + damfc;
654    }
655    load = HystereticLoad::Negative;
656
657    if rot_min < HYSTERETIC_NEG_INF_STRAIN {
658        rot_min = HYSTERETIC_NEG_INF_STRAIN;
659    }
660    rot_min = if rot_min < rot1n { rot_min } else { rot1n };
661
662    let minmom = envelope.neg_stress(rot_min);
663    let rotlim = envelope.pos_rotlim(rot_max);
664    let rotrel = if rotlim < rot_pu { rotlim } else { rot_pu };
665
666    let rotmp2 = rot_min - (1.0 - pinch_y) * minmom / (eun * kn);
667    let rotch = rotrel + (rotmp2 - rotrel) * pinch_x;
668
669    let (stress, tangent) = if strain > rot_pu {
670        let t = eup * kp;
671        let s = cstress + t * dstrain;
672        if s <= 0.0 {
673            (0.0, eup * 1.0e-9)
674        } else {
675            (s, t)
676        }
677    } else if strain > rotch {
678        if strain >= rotrel {
679            (0.0, eun * 1.0e-9)
680        } else {
681            let t = minmom * pinch_y / (rotch - rotrel);
682            let tmpmo1 = cstress + eun * kn * dstrain;
683            let tmpmo2 = (strain - rotrel) * t;
684            if tmpmo1 > tmpmo2 {
685                (tmpmo1, eun * kn)
686            } else {
687                (tmpmo2, t)
688            }
689        }
690    } else {
691        let t = (1.0 - pinch_y) * minmom / (rot_min - rotch);
692        let tmpmo1 = cstress + eun * kn * dstrain;
693        let tmpmo2 = pinch_y * minmom + (strain - rotch) * t;
694        if tmpmo1 > tmpmo2 {
695            (tmpmo1, eun * kn)
696        } else {
697            (tmpmo2, t)
698        }
699    };
700
701    (stress, tangent, rot_min, rot_pu, load)
702}
703
704#[cfg(test)]
705mod tests {
706    use super::*;
707
708    /// A symmetric, 3-point-per-side, bilinear-then-flat backbone shared by
709    /// every `Hysteretic` test below: `mom1p/rot1p=10/0.01` (E1p=1000),
710    /// `mom2p/rot2p=15/0.02` (E2p=500), `mom3p/rot3p=15/0.03` (E3p=0, flat
711    /// plateau) — negative side the exact mirror. `beta=0` makes the
712    /// ductility-degradation factors `kn`/`kp` identically `1.0`
713    /// throughout (`x^0 == 1`, and `1.0 < 1.0` is false so the `1.0/kn`
714    /// branch never engages either) — deliberately excluded from these
715    /// tests to keep the hand computation tractable; `damfc1=damfc2=0`
716    /// disables the deformation/energy damage factor for the same reason.
717    fn hysteretic_test_material() -> Material {
718        Material::hysteretic(
719            10.0, 0.01, 15.0, 0.02, 15.0, 0.03, -10.0, -0.01, -15.0, -0.02, -15.0, -0.03, 0.5, 0.5,
720            0.0, 0.0, 0.0,
721        )
722    }
723
724    #[test]
725    fn hysteretic_gives_exact_initial_tangent_at_zero_strain() {
726        let m = hysteretic_test_material();
727        assert_eq!(m.trial_stress_tangent(0.0), (0.0, 1000.0));
728    }
729
730    #[test]
731    fn hysteretic_matches_bilinear_envelope_on_first_loading() {
732        let m = hysteretic_test_material();
733        // Within the first segment (strain <= rot1p=0.01): E1p=1000.
734        assert_eq!(m.trial_stress_tangent(0.005), (5.0, 1000.0));
735        // Within the second segment (rot1p < strain <= rot2p=0.02):
736        // mom1p + E2p*(strain-rot1p) = 10 + 500*0.005 = 12.5.
737        assert_eq!(m.trial_stress_tangent(0.015), (12.5, 500.0));
738    }
739
740    #[test]
741    fn hysteretic_pinches_on_reload_after_a_full_reversal() {
742        let m = hysteretic_test_material();
743        // Push to the positive envelope's second breakpoint (rot_max
744        // becomes exactly rot2p=0.02, stress=mom2p=15), then reverse all
745        // the way past the negative breakpoint too (rot_min becomes
746        // exactly rot2n=-0.02, stress=mom2n=-15) — both reach the
747        // envelope-extreme branch directly, since rot_min/rot_max both
748        // start at 0.0 (any first excursion in either direction is "new").
749        let m = m.commit(0.02).commit(-0.02);
750
751        // Now reload from -0.02 toward 0.02 — strictly interior to
752        // [rot_min, rot_max], so this hits `hysteretic_positive_increment`
753        // for the first time, which captures the pinching reference point
754        // `rot_nu = cstrain - cstress/(Eun*kn) = -0.02 - (-15)/1000 =
755        // -0.005` off the just-completed reversal.
756        //
757        // At strain=-0.01 (still below rot_nu=-0.005): pure elastic
758        // reloading at slope Eun*kn=1000 from the reversal point:
759        // stress = cstress + 1000*dstrain = -15 + 1000*0.01 = -5.
760        let (stress, tangent) = m.trial_stress_tangent(-0.01);
761        assert!((stress - -5.0).abs() < 1e-9, "got {stress}");
762        assert!((tangent - 1000.0).abs() < 1e-9, "got {tangent}");
763
764        // At strain=0.0 (past rot_nu, in the pinched-plateau region below
765        // rotch): the pinched-target line `(strain-rotrel)*maxmom*pinchY/
766        // (rotch-rotrel)` governs rather than the steep elastic-unload
767        // line — the hallmark of pinching (a much softer path back through
768        // zero than straight elastic reload would give).
769        let (stress, tangent) = m.trial_stress_tangent(0.0);
770        assert!((stress - 4.285714285714286).abs() < 1e-9, "got {stress}");
771        assert!((tangent - 857.1428571428571).abs() < 1e-9, "got {tangent}");
772        assert!(
773            tangent < 1000.0,
774            "pinched stiffness should be softer than the elastic Eun, got {tangent}"
775        );
776    }
777
778    #[test]
779    fn hysteretic_energy_accumulates_as_half_trapezoid_of_stress_over_strain() {
780        let m = hysteretic_test_material();
781        let m = m.commit(0.01); // exactly rot1p: stress=mom1p=10.
782        let Material::Hysteretic(fields) = m else {
783            panic!()
784        };
785        let energy_d = fields.energy_d;
786        // Triangular area under the elastic segment: 0.5*10*0.01 = 0.05.
787        assert!((energy_d - 0.05).abs() < 1e-9, "got {energy_d}");
788    }
789}
carapace_core/model/materials/hysteretic.rs

carapace_core/model/materials/
hysteretic.rs
