Add Modular UMAT integration#61
Open
chemiskyy wants to merge 53 commits into
Open
Conversation
Include 'simcoon.modular' in the conda.recipe/meta.yaml test.imports list so the conda build/test verifies the modular subpackage can be imported. This ensures the modular component is packaged and available during CI or conda package validation.
Introduce a new DamageMechanism implementing isotropic scalar continuum damage mechanics. Adds header and implementation files that integrate with StrainMechanism and InternalVariableCollection, register history variables (D, Y_max), compute the energy-based driving force Y = 0.5 * sigma : S : sigma, and evaluate damage via linear, exponential, power-law, and Weibull evolution laws. Includes damaged stiffness/tangent contributions, Jacobian/work computations, parameter configuration and basic validation, and caching of the compliance tensor.
Introduce PlasticityMechanism class (header + implementation) to provide a rate-independent plasticity strain mechanism for the modular UMAT. The new files combine a yield criterion, isotropic hardening and kinematic hardening (via factory-created hardening objects), register internal variables (p, EP and backstresses), and implement core routines: configure, register_variables, compute_constraints, compute_flow_directions, compute_jacobian_contribution, inelastic_strain, update, tangent_contribution and compute_work. The implementation caches flow direction and related quantities for Jacobian/tangent calculations, supports multiple yield/hardening types, and uses Armadillo vectors/matrices and the existing InternalVariableCollection API.
Introduce ElasticityModule (header + implementation) to encapsulate linear elastic behavior for the modular UMAT. Supports isotropic, cubic (EnuG or Cii conventions), transversely isotropic, and orthotropic configurations via configure_* helpers and a configure(type, props, offset) convenience. Stores 6x6 stiffness (L) and compliance (M) matrices, 6-component CTE vector (alpha), and provides derived helpers damaged_L(d) and thermal_strain(DT). Uses Armadillo types and existing constitutive functions (L_iso/M_iso, L_cubic/M_cubic, L_isotrans/M_isotrans, L_ortho/M_ortho). Adds basic validation (configured_ flag) and a constexpr props_count helper.
Introduce modular hardening framework: new header and implementation for isotropic and kinematic hardening models. Implements No/Linear/Power-law/Voce/Combined-Voce isotropic models and No/Prager/Armstrong–Frederick/Chaboche kinematic models, with factory create() functions, property parsing (configure), evaluation (R, dR/dp), internal-variable registration, backstress accumulation, flow rules and update routines. Uses Armadillo and InternalVariableCollection for vector storage and exposes props_count helpers to consume material property arrays.
Add a type-safe InternalVariable class to support constitutive-model state storage. New files include/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable.hpp and src/Continuum_mechanics/Umat/Modular/internal_variable.cpp implement scalar, 6-component Voigt vector, and 6x6 matrix types with start/current state tracking, delta computation, rotation/objectivity (uses rotate_strain / rotateL), and pack/unpack serialization to a flat statev array with configurable offset. Uses Armadillo, performs input validation, and throws on type mismatches. Intended for use in UMAT modular constitutive implementations.
Introduce InternalVariableCollection (header and implementation) to manage multiple InternalVariable instances for modular UMAT models. Provides named registration (scalar/6-vector/6x6-matrix), automatic statev offset computation, pack/unpack, rotation, start/reset helpers, iteration, and utility functions (names, clear). Uses unique_ptr storage (non-copyable, moveable), throws on invalid access/duplicate names, and depends on Armadillo.
Introduce a composable ModularUMAT class and standalone umat_modular function to orchestrate elasticity and multiple strain mechanisms (plasticity, viscoelasticity, damage). Adds public configuration APIs (set_elasticity, add_plasticity/add_viscoelasticity/add_damage, configure_from_props), initialization and state handling (internal variable collection), the main run() entry that performs elastic prediction, return mapping (Newton + Fischer–Burmeister), consistent tangent computation, and work calculations. Implements compute_tangent and return_mapping, and integrates with existing mechanism classes and numeric solver utilities. Files added: include/.../modular_umat.hpp and src/.../modular_umat.cpp (GPL header included).
Introduce a new header defining the StrainMechanism abstract base class for modular UMATs. Adds MechanismType enum and a comprehensive interface for pluggable strain mechanisms (configuration, variable registration, constraint evaluation, flow directions, Jacobian contributions, inelastic strain, internal variable updates, tangent contributions, and work decomposition). Uses Armadillo and InternalVariableCollection and includes GPL license header; intended to standardize integration of plasticity, viscoelasticity, damage, etc.
Introduce ViscoelasticMechanism implementing a generalized Maxwell (Prony series) viscoelastic branch model. Adds header and source with configuration, internal-variable registration, constraint evaluation, flow directions, Jacobian contribution, inelastic strain accumulation, update routine, tangent modification and work computation. Supports N Prony terms, registers one viscous strain vector per branch, reads g_i and tau_i from property vector, and uses an optional reference stiffness (L_0) / compliance for branch targets. Includes basic runtime checks (positive tau) and a simple algorithmic tangent and update/constraint approximations (implicit exponential integration and norm-based constraints). Files use Armadillo and include GPL license headers.
Introduce YieldCriterion class (header and implementation) to provide a modular wrapper around existing Eq_stress/dEq_stress criterion functions for UMAT. Adds YieldType enum (VON_MISES, TRESCA, DRUCKER, HILL, DFA, ANISOTROPIC), configuration helpers (including configure from a props vector with offset), parameter storage via arma::vec, and public APIs to compute equivalent_stress, flow_direction and plastic_flow (with optional backstress shift). Includes a constexpr props_count utility and runtime checks that the criterion is configured; integration relies on simcoon/Continuum_mechanics/Functions/criteria.hpp.
Introduce a new simcoon.modular module implementing a composable UMAT material system. Adds enums and dataclasses for elasticity (isotropic, cubic, transverse isotropic, orthotropic), yield criteria (VonMises, Tresca, Drucker, Hill, DFA, anisotropic), isotropic and kinematic hardening laws (linear, power-law, Voce, combined Voce, Prager, Armstrong-Frederick, Chaboche), and mechanisms (Plasticity, Viscoelasticity with Prony series, Damage). Implements ModularMaterial to assemble a flat props array (for the C++ "MODUL" UMAT), compute nprops/nstatev, provide human-readable summaries, and includes convenience factory functions (elastic_model, elastoplastic_model, viscoelastic_model).
Add an explicit `from simcoon import modular` to simcoon.__init__ so the `simcoon.modular` submodule is available to consumers when importing the package. This makes the modular API accessible without requiring an extra import.
Delete generated test output files from testBin/Umats/UMEXT/results: results_job_global-0.txt and results_job_local-0.txt. These are result artifacts (deleted mode 100755) and are not intended to be tracked in source control.
Add comprehensive GTest unit tests for the Modular UMAT implementation (test/Libraries/Umat/Modular/Tmodular_umat.cpp) and register the new test in test/CMakeLists.txt. The tests cover InternalVariable and InternalVariableCollection, ElasticityModule (isotropic/cubic/CTE), YieldCriterion (von Mises), Isotropic and Kinematic hardening models, PlasticityMechanism behavior, ModularUMAT initialization/run behavior, and an integration test for uniaxial tension. This provides broad verification of elastic/plastic responses, packing/unpacking of state variables, hardening laws, and consistency of tangents and stresses.
Integrate the new modular UMAT into the project: add modular_umat includes to the Python wrapper and smart UMAT selector, register the "MODUL" UMAT id in both mappings (30 for Python wrapper, 200 for smart selector), and add case branches to call umat_modular accordingly. Update CMakeLists to include the Modular UMAT source files (internal_variable, elasticity_module, yield_criterion, hardening, plasticity_mechanism, viscoelastic_mechanism, damage_mechanism, and modular_umat) so the modular implementation is built. This wires the composable constitutive model into the build, runtime selection, and Python bindings.
Brings the modular UMAT framework (ElasticityModule, PlasticityMechanism, ViscoelasticMechanism, DamageMechanism, hardening, yield criteria, InternalVariable/Collection, ModularUMAT orchestrator, MODUL UMAT name) on top of the Tensor2/Tensor4/Rotation/Fastor stack. Conflict resolutions: - conda.recipe/meta.yaml: keep both simcoon._core and simcoon.modular imports - umat_smart.cpp / Python umat.cpp: merge SMA name expansion + MODUL registration - internal_variable.cpp: replace removed rotateL() with Rotation::from_matrix(DR).apply_stiffness(L) - *_mechanism.cpp / hardening.cpp / damage_mechanism.cpp: rename sim_iota to simcoon::iota (renamed during the Tensor migration) Baseline tests pass: 63 Ttensor + 27 Tmodular_umat.
Adds typed Tensor2/Tensor4 accessors alongside the raw arma::vec/arma::mat accessors so callers can use the type-safe Tensor API (mirrors the EPICP plugin pattern): - L_tensor() -> Tensor4 (Tensor4Type::stiffness) - M_tensor() -> Tensor4 (Tensor4Type::compliance) - alpha_tensor() -> Tensor2 (VoigtType::strain) - damaged_L_tensor(d) -> Tensor4 (stiffness) - thermal_strain_tensor(DT) -> Tensor2 (strain) Backward-compatible: L(), M(), alpha(), damaged_L(), thermal_strain() unchanged. Internal storage stays arma::mat/vec; the typed accessors construct on demand. Adds ElasticityModuleTest.TensorAccessors verifying numerical equivalence and correct type tags, plus an end-to-end stress = L : eps round-trip through the typed contract.
Adds typed views over the stored Voigt vector / 6x6 matrix: - as_tensor2(VoigtType) / as_strain() / as_stress() - as_tensor4(Tensor4Type) / as_stiffness() / as_compliance() - set_tensor2(t) / set_tensor4(t) for write-back The caller picks the VoigtType matching what the stored variable represents (plastic strain → strain, etc.) since the storage itself is untyped. The default VoigtType::strain matches the existing rotate() convention which already uses the strain-Voigt factor-2 layout. Also adds rotate(const Rotation&) overload so callers with a Rotation object don't need to round-trip through a 3x3 matrix. Tests: TensorViews (round-trip + type-mismatch rejection) and RotationOverload (Rotation overload matches mat overload bitwise).
Adds typed overloads: - equivalent_stress(const tensor2&) -> double - equivalent_stress(const tensor2&, const tensor2&) -> double (with backstress) - flow_direction(const tensor2&) -> tensor2 (VoigtType::strain) - flow_direction(const tensor2&, const tensor2&) -> tensor2 (with backstress) Implementations delegate to the existing arma::vec overloads via .voigt(). The flow direction returns strain-Voigt (factor-2 on shear) to match the simcoon convention where dEP = dp * Lambda updates plastic strain additively in the strain-Voigt storage layout. Test: VonMisesTensor2Match verifies bitwise equivalence with the arma overloads and the correct VoigtType on the returned flow direction.
…back-strain refactor C++ side of the modular UMAT framework: - ViscoelasticMechanism: faithful port of Prony_Nfast kernel — per-branch L_i(E_i, nu_i) and H_i(etaB_i, etaS_i), strain-form equality constraint Phi_i = ||invH_i.L_i.(eps - EV_i)|| - dv_i/dt = 0, lead scalar v_i + EV_i per branch (7 statev each). Inelastic strain via sum_i (M_0.L_i).EV_i. - StrainMechanism: cross-mechanism coupling via three new virtuals: dPhi_dsigma(), kappa(), K_cross(). Rate-dependent mechanisms with strain- form Phi return empty dPhi_dsigma; FB residual handles equality+complement uniformly. - ModularUMAT::return_mapping reorganised into three phases: (1) per-mechanism constraint evaluation, (2) global B_lj = -dPhi^l/dsigma . kappa^j + K^lj cross assembly, (3) per-mechanism diagonal fill. - IVC prefix-based namespacing (plast<i>_, visco<i>_, damage<i>_) lets multiple same-type mechanisms coexist; StrainMechanism::set_ivc_prefix propagates to KinematicHardening for the alpha keys. - Backstress refactored: store back-strain alpha (true thermodynamic internal variable, strain-like Voigt) instead of X (derived force). total_backstress() returns X = (2/3) C alpha derived. Matches the original plastic_chaboche_ccp.cpp convention. file-local helper backstress_t(C, alpha) centralises the conjugacy. - InternalVariable: vtype_ (was rotation_kind_) is single VoigtType source of truth; rotate(Rotation&) routes through tensor2.rotate(R) so strain/stress dispatch is in one place. as_tensor2() defaults to stored vtype_. - Mechanism kernels in hardening / damage migrated to typed tensor2 / tensor4 algebra where natural; viscoelastic stays raw arma because L_i / H_i are not modelled by tensor4. - Hot-path: cached IVC keys at register_variables time across all mechanisms; cached mech_offset_ on ModularUMAT; dPhi_dsigma/kappa return const-ref into mechanism-local cache buffers; eliminates ~5-8x of allocator pressure in the FB hot loop. - Tests (Tmodular_umat.cpp): TwoPlasticityEquivalentToOne, TwoPlasticityMechanismsInitialize, ArbitraryMultiMechanismInitialize, ViscoelasticRelaxation, ChabocheDuplicateEqualsDoubledAF, CombinedVoceDuplicateEqualsDoubled, RotationOverload, TensorViews, YieldAccessors. Total: 38 tests passing.
- simcoon/modular.py:
* Viscoelasticity dataclass now takes (E, nu, etaB, etaS) per Prony
branch (was (g, tau)) — matches the Prony_Nfast C++ kernel.
nstatev = 7 * N (was 6 * N): adds the lead scalar v_i per branch.
* __post_init__ rejects the legacy 2-tuple form and AF/Prager called
with list arguments, with pointed error messages.
* ModularMaterial.summary() prints per-branch viscoelastic parameters
instead of just "N Prony terms".
* viscoelastic_model() factory updated to the new 4-tuple layout.
- examples/umats/MODUL.py: declarative-API end-to-end example
(ModularMaterial -> sim.solver -> stress-strain plot). Uses isotropic
elasticity + von Mises + Voce hardening on a monotonic ramp + 2 cyclic
blocks.
- examples/umats/data/MODUL_path.txt: companion path file.
- test_modular.py: pytest regression for the modular Python path.
test_two_plasticity_equivalent_to_one_linear protects the FB cross-
Jacobian fix (would have caught the pre-fix 19000 MPa nonphysical
stress regression). Three other tests cover input validation and
Voce smoke.
Conflict resolutions:
- pybind umat.cpp: MODUL renumbered 30 -> 200 (feature/Tensor took 30 for
NEOHI); map + dispatch case updated to match the C++ id
- umat_smart.cpp: feature/Tensor map (SMRDI..=9, EPTRI=17) + {MODUL,200}
- umat_modular gains trailing tangent_mode param (unused) to match the
unified small-strain UMAT pointer signature
- test/CMakeLists.txt: union (Ttangent_* + Tmodular_umat)
- MODUL example relocated to examples/mechanical + examples/data, test
paths updated accordingly
tensor4 now stores Kelvin-Mandel internally; mat() recomputes the engineering form, exact only to ~1 ulp of the entries (~2.5e-11 abs for a 210 GPa stiffness). 1e-12 was calibrated to the old zero-conversion storage; 1e-9 keeps the check meaningful (1e-14 relative).
tangent_contribution rebuilt K(i,i) as -dot(dPhi_i_dv, kappa_i) alone, dropping the -1/DTime term that compute_constraints caches in K_diag_. The Prony_Nfast reference keeps that term (Bhat = -K); without it the rank-one correction roughly doubles and Lt turns singular/indefinite, so any mixed-control solver run diverged (inv(): matrix is singular). Reuse the cached K_diag_ instead. Add test_viscoelastic_matches_pronk_reference: MODUL 3-branch Prony vs the reference PRONK UMAT through sim.solver on PRONK_path.txt — max stress deviation now 1.8e-4 relative (was: divergence). Closes the missing viscoelastic end-to-end coverage.
StrainMechanism::dPhi_dsigma()/kappa() now return const std::vector<tensor2>& so each cross-Jacobian operand carries its Voigt convention: dPhi strain-typed (dEq_stress), kappa stress-typed for plasticity/viscoelasticity, and strain-typed for damage (M.sigma product) — the damage convention mix in the B assembly is now visible in the type system and documented as a deliberate modeling item. The orchestrator contracts on engineering Voigt components (dot of .voigt()), the work-conjugate pairing — bit-identical to the previous raw dots (strain voigt round-trips /2 then *2 exactly in FP). Also hoist the per-mechanism kappa() fetch out of the phase-2 row loop (it was re-invoked for every constraint row). Goldens (elastoplastic/viscoelastic/combined, 600 incr each): outputs bit-identical, timings within noise. gtests 38/38, pytest 5/5.
L_tensor()/M_tensor() built a tensor4 (eng->Mandel congruence) on every call; with the typed plasticity core about to contract L_tensor() inside the FB loop that would be a per-iteration cost. Build the typed mirrors once in refresh_tensors(), called by every configure_* path; accessors now return const&. Goldens bit-identical, timings in noise.
…hmarked) Add KinematicHardening::total_backstress_t() (strain-convention tensor2, documenting the back-strain Voigt layout the raw accessor returns) and tensor2 convenience overloads of hardening_modulus()/update() that delegate to the raw virtuals. Typing the PlasticityMechanism hot members (flow_dir_/kappa_ as tensor2, EPICP-style) was tried and REVERTED: the FB return-map re-derives the flow every iteration against raw yield/hardening internals, so typed members forced eng<->3x3 conversions at each boundary — a reproducible +7% on the elastoplastic golden (vs noise +-3%). The typed view is produced once per iteration at the dPhi_dsigma()/kappa() interface instead; rationale recorded in plasticity_mechanism.hpp. This is the 'typed cores, raw solver' boundary the EPICP plugin pattern implies for a per-iteration CCP loop. Goldens bit-identical; timings back in noise. gtests 38/38, pytest 5/5.
YieldCriterion tensor2 overloads no longer round-trip through .to_arma_voigt() for VON_MISES: equivalent_stress uses Mises(tensor2) and flow_direction uses flow_normal(tensor2), both evaluated directly on the 3x3 (the kernels added for the EPICP typed plugin). Parametric criteria (Drucker = J2-J3 form, not Drucker-Prager; Hill/DFA/Ani) keep delegating to the raw Eq_stress dispatch. Tests: DruckerTensor2Delegates pins the delegate equivalence; VonMisesNearZeroStress pins the cone-vertex behavior and documents the intentional difference from the raw path there (flow_normal zeroes the normal below iota, eta_stress normalizes noise for any n > 0). gtests 40/40, pytest 5/5, goldens bit-identical.
…tests DamageMechanism evolved D and softened the TANGENT ((1-D)*Lt in tangent_contribution) but the orchestrator's stress prediction ignored damage entirely — el_pred always used the undamaged L, so the composed response never softened and Newton saw an inconsistent tangent/residual pair. Add StrainMechanism::stiffness_reduction() (default 1; damage returns 1-D; factors compose multiplicatively) and scale the elastic prediction sigma = f * L : Eel at both return_mapping stress sites. New coverage for the previously untested composed-damage paths: - gtest DamageElasticUniaxial: D monotone in [0, D_c], stress matches (1-D)*L:E within 5%, D frozen on unloading (Y_max history) - gtest PlasticityDamageCombined: plasticity+damage both activate (exercises phase-2 cross-assembly with the strain-typed damage kappa), softened vs plasticity-only reference - pytest test_damage_softens_end_to_end through sim.solver gtests 42/42, pytest 6/6 (test_core 268), goldens bit-identical, timings in noise.
Move conceptual documentation into headers, per house style (.hpp = Doxygen contract, .cpp = load-bearing implementation notes only): - DamageType enum now carries the four evolution-law formulas in LaTeX (with the Y_eff = max(Y, Y_max) / threshold / D_c cap semantics); drop the .cpp per-case restatements - PragerHardening class doc gains the back-strain convention block (alpha is the stored internal variable, X = (2/3) C alpha is derived) formerly an implementation comment in hardening.cpp - Delete narration one-liners in elasticity_module.cpp (Compute stiffness/compliance/CTE) and the factory-method banners Behavior-neutral: gtests 42/42, pytest 6/6.
…anch
conda.recipe/meta.yaml selects its source via CONDA_BUILD_LOCAL:
git_url: ... # [not environ.get('CONDA_BUILD_LOCAL')]
path: .. # [environ.get('CONDA_BUILD_LOCAL')]
No workflow ever set it, so every conda-packaging run cloned the
repository default branch (master) and built THAT against the PR's
recipe. PR #61's recipe adds 'import simcoon.modular' to the test
imports -> all three OS jobs failed with ModuleNotFoundError even
though the module exists on the branch (the built package even shipped
master-only FTensor headers, confirming the wrong source).
Set CONDA_BUILD_LOCAL=1 at the job level: PR runs now validate the PR
source, and the workflow_call 'ref' input actually affects what is
built.
tangent_mode now reaches the orchestrator (umat_modular -> run -> compute_tangent; umat_smart passes umat_M->tangent_mode). Mode 1 assembles the Simo-Hughes consistent tangent of the coupled sub-system through master's multi-mechanism assemble_algorithmic_tangent(): the phase-2/3 Jacobian assembly is factored into assemble_jacobian(), the opted-in sub-block is handed over as Bhat = -B with the mechanisms' kappa / dPhi_dsigma caches and the new flow Hessians. Mechanisms whose flow is stress-independent (Prony viscoelastic, scalar damage) keep their continuum tangent_contribution in every mode - exact, not a fallback. New interface (also the tangent_mode==2 preparation, aligned with return_mapping.hpp on feature/algorithmic_tangent): - StrainMechanism::dLambda_dsigma(): flow Hessians per constraint, compliance-typed tensor4 (stress->strain map: Mandel congruence carries the Hessian shear factors by construction). Plasticity implements it at the shifted stress via YieldCriterion::flow_hessian (deta_stress / ddHill_stress / ddDFA_stress / ddAni_stress); Tresca and Drucker have no analytic Hessian and opt out. - StrainMechanism::refresh_state(): backward-Euler state rebuild from the IVC start values under total Ds (the CPP update_state contract). Plasticity implements it with the AF/Chaboche closed forms (alpha = (alpha_n + dp n)/(1 + D dp)) and a fixed point on the n <-> X coupling; KinematicHardening::refresh_state added. - ModularUMAT::set_solver_params(maxiter, precision). BUG FIX exposed by the mode-1 finite-difference test: the plasticity local Jacobian carried +H_total instead of -H_total (dPhi/dp = -H; the plastic_isotropic_ccp convention is K = -H, B = -n:kappa - H). The wrong Newton slope skewed the FB error metric, so the loop reported convergence prematurely: 4.4% peak-stress error vs the EPCHA reference on the cyclic Voce+Chaboche path (31.9 MPa at 728 MPa peak), invisible to all monotonic tests. Fixed, MODUL now matches EPCHA to 3.7e-4 relative over 600 cyclic increments. Tests: AlgorithmicMatchesFiniteDifference (mode-1 tangent = FD Jacobian of the discrete update to <1e-3, continuum >10x worse), RefreshStateBackwardEulerAF (implicit relation satisfied to 1e-10), test_tangent_mode_1_same_converged_response (solver mode 0 == mode 1), test_chaboche_matches_epcha_reference (regression for the sign bug). gtests 44/44, pytest 8/8 (test_core 273). Goldens re-recorded post-fix (old elastoplastic golden was 4.4% off the reference; combined case now converges 3x faster with the honest slope).
The header documented a variable-term layout (props[4]=N_iso, props[5]=N_kin, per-term tables, nstatev = 8 + 6 N_kin = 20) that the code never implemented: umat_plasticity_chaboche_CCP reads a FIXED 10-prop layout (E, nu, alpha, sigmaY, Voce Q, b, C_1, D_1, C_2, D_2 — one Voce term in rate form dHp = b(Q-Hp) dp, exactly two AF backstresses) and 33 statev (T_init, p, EP, back-strains a_1/a_2, backstresses X_1/X_2, Hp). Following the stale table (nstatev=20) aborts with an out-of-bounds statev read. Rewritten: props table (fixed indices), statev table (full 33-entry layout, a_i vs X_i distinction), H_tot formula, Voce rate form, and the @code example (props(10), zeros(33)). Points users needing a variable number of terms to the modular UMAT (MODUL).
…n; IV hardening BUG FIX (found by the review's typed-rotation test, which built its fixture by copying the module's call): configure_orthotropic passed (E1, nu12, nu13, E2, nu23, E3, ...) to L_ortho/M_ortho, whose EnuG slot order is (E1, E2, E3, nu12, nu13, nu23, ...) — Poisson ratios landed in the Young's-moduli slots, producing an indefinite stiffness (cond ~ 3.5e9, negative eigenvalue) for every orthotropic modular material. The props chain (Python dataclass -> configure -> named params) was correct; only the call site was wrong. Orthotropic was the one elasticity type without a configuration test — added OrthotropicConfiguration (match vs direct L_ortho + positive definiteness). InternalVariable review batch: - MATRIX_6x6 rotation now routes through tensor4 with a new stored Tensor4Type (t4type_, ctor/add_mat param, default stiffness) — the correct congruence for stiffness- AND compliance-like storage (Matrix66TypedRotation pins M_rot == inv(L_rot) for orthotropic L); as_tensor4() defaults to the stored type, mirroring as_tensor2(). - set_tensor2/set_tensor4 re-express the value in the variable's own convention (via the 3x3 / Mandel form) instead of copying raw components — assigning a stress-typed tensor2 to a strain-typed variable can no longer silently drop the factor-2 shear convention (SetTensor2ReexpressesConvention). - pack/unpack throw on out-of-bounds offsets instead of silently truncating (PackUnpackOutOfBoundsThrows). - rotate_all builds the Rotation once instead of re-extracting the quaternion from the same 3x3 for every variable. gtests 48/48, pytest 273, goldens bit-identical.
ElasticityModule review batch: - refresh_tensors() now rejects a non-positive-definite stiffness via a Cholesky test (cheap SPD check; umat_modular rebuilds the module every integration-point call, so a full eig_sym there was measurable). This turns an inadmissible material - a wrong parameter order (cf. the orthotropic arg-order bug) or physically impossible constants (e.g. nuTL^2 * EL/ET too large for transverse isotropy) - into an immediate, attributable error instead of an indefinite L that fails deep in the return mapping. - configure_transverse_isotropic validates the axis (1/2/3) BEFORE calling L_isotrans, which otherwise terminates the process with exit(0) on an invalid axis - uncatchable from Python or a solver. - Doc note: the cubic Cii form is reachable only through the direct C++ API; the props-driven configure() and the Python dataclasses use EnuG. Tests: RejectsInadmissibleMaterial (indefinite ortho, invalid axis, inadmissible transverse-iso all throw; admissible carbon-fibre-like constants accepted). gtests 49/49, pytest 273, goldens in noise. Note: damaged_L / damaged_L_tensor are now exercised only by their own unit test - the damage mechanism computes (1-D)L inline and the orchestrator uses stiffness_reduction(). Kept as public convenience accessors.
yield_criterion review: the criterion params (Drucker b/n, Hill FGHLMN, DFA +K, Ani 9-term) were verified against criteria.cpp P_Hill/P_DFA/ P_Ani layouts (all correct), but no anisotropic criterion had ANY end-to-end modular coverage — the same gap that hid the orthotropic arg-order bug. Add test_hill_matches_ephil_reference: MODUL Hill+linear hardening == EPHIL to machine precision (max |dsig11| = 0). Uses linear hardening (m=1) on purpose: the power-law tangent dR/dp = m k p^(m-1) is singular at p=0, so with m<1 the two CCP loops differ by a re-converging transient (~1.4% at yield onset, verified harmless); m=1 removes it and isolates the Hill criterion, which is exact. pytest 9/9.
…rror) BUG #6 (pre-existing, found by the review's shear cross-check vs EPCHA): the backstress X = (2/3) C a is a STRESS-like tensor - it shifts the stress in the yield function (sigma - X) and pairs with the strain-like flow n in X:n - but it was carried in the back-strain (engineering strain) Voigt convention with a factor-2 on shear. Under shear the residual therefore subtracted a doubled X: 24.7% peak-stress error vs the reference EPCHA (sig12 287 vs 230 MPa). Invisible on the uniaxial path (shear components zero), so the uniaxial EPCHA test passed at 3.7e-4 and hid it - and isotropic von Mises under shear was already bit-identical, pinning the defect to the kinematic term. Two sites, both needed: - backstress_t(): re-tag (2/3) C a as stress (via the convention-free 3x3) so its Voigt form has no factor-2 shear. - ChabocheHardening::total_backstress(): accumulate STRESS-typed; a strain-typed zero accumulator re-tagged the sum and put the factor-2 back (this is the single-AF-via-Chaboche path too). hardening_modulus() X:n and compute_work() X:dEP are fixed for free (both now dot a stress-typed X against strain-typed n / dEP). Now MODUL Chaboche under shear matches EPCHA to 3.5e-5. New regression test_chaboche_shear_matches_epcha_reference on a pure-shear path (E12 driven, 2 backstresses) + examples/data/SHEAR_path.txt. Uniaxial goldens bit-identical (zero shear). gtests 49/49, pytest 10/10.
…flow_directions) strain_mechanism review — two dead-API removals, zero callers anywhere (orchestrator, tests, pybind, Python): - active_ / is_active() / set_active(): a settable flag the orchestrator never reads. set_active(false) silently left the mechanism fully contributing — worse than dead code (a misleading no-op). Toggling a mechanism needs real offset/statev handling; removed rather than left as a trap. - compute_flow_directions(): pure virtual + 3 implementations, never invoked. Flow directions are computed inline in compute_constraints (cached as flow_dir_ / Lambda_i_); this was vestigial from an earlier design where they were a separate pass. Behavior-neutral. gtests 49/49, pytest 10/10.
… vars + docs viscoelastic_mechanism: correct the misleading 'filled in a separate commit' note — the constraint Jacobian is diagonal BY STRUCTURE (Phi_i depends on EV_i alone; dPhi_i/dv_j = 0 for j != i). Branch coupling enters only the consistent tangent. Validated by the PRONK equivalence test (3-branch, 1.8e-4). damage_mechanism: - remove dead D_new (computed, never used) and the unused dY multiplier read in update(). - document what the code actually does: damage is integrated EXPLICITLY (update() sets Y_max = max(Y, Y_max), D = f(Y_max)); the FB multiplier is unused and B(r,r) = 1.0 is a nominal unit slope for the history-type row, which self-satisfies (Phi = Y - Y_max -> 0 after the Y_max update). Cross-mechanism plasticity<->damage coupling still flows through the orchestrator's off-diagonal dPhi_dsigma . kappa terms. Behavior-neutral. gtests 49/49, pytest 10/10.
The orchestrator committed whatever the FB loop produced, even a diverged NaN/Inf state, without signalling the global solver — a hard increment would silently propagate garbage (and never retry). Now run() checks sigma.is_finite() after return_mapping and, on failure, sets tnew_dt = 0.5 (halve the increment and retry), skips pack_all so statev stays at its incoming values for the retry, and sets Lt elastic to keep the rejected-step output well-defined. Trigger is non-finite, NOT error > precision_: a FINITE unconverged-at- maxiter state is still committed, matching the reference CCP UMATs (UNUSED(tnew_dt)) — in particular the damage row is integrated explicitly and never drives the FB error to precision_, so an error-based cut would wrongly reject every damage step. Behavior-neutral on all converging cases: gtests 49/49, pytest 275, goldens bit-identical.
required_nstatev() hand-counts the statev size per mechanism type so callers can size their statev array before initialize(). That duplicates the authoritative count register_variables() produces, and could drift silently (add a variable to a mechanism, forget the hand-count -> a too-small statev and corrupted state downstream). initialize() already has both numbers; it now throws if they disagree, so any drift fails loudly on the first modular run instead of much later. gtests 49/49, pytest 10/10.
…der shear modular.py review: props/nstatev assembly matches configure_from_props / register_variables for every dataclass (verified — and now guarded at runtime by the initialize() consistency check). The one kinematic path not covered by an equivalence test was the ArmstrongFrederickHardening dataclass, which routes to a DISTINCT C++ class (not Chaboche). Lock it: AF == single-term Chaboche bit-for-bit under shear, exercising the AF total_backstress path (backstress_t, no accumulator). pytest 11/11.
Add runtime validation for elasticity configuration and corresponding tests. Include <stdexcept> and check that the assembled stiffness matrix L_ is strictly positive definite in refresh_tensors(), throwing a runtime_error if not (to prevent downstream failures from indefinite stiffness). Validate the transverse-isotropic axis in configure_transverse_isotropic() and throw on invalid axis instead of allowing L_isotrans to terminate the process. Minor refactor of a switch default comment and a note about cubic Cii API surface. Update unit tests to assert that inadmissible material configurations and invalid axis values are rejected (EXPECT_THROW) and that valid axis usage does not throw.
Remove deprecated damaged-stiffness helpers and update tests. Deleted the damaged_stiffness static method from DamageMechanism (header + source) and removed damaged_L / damaged_L_tensor from ElasticityModule (header + source). Updated unit test to drop assertions that relied on the removed API. This is an API cleanup — no other behavior changes are introduced.
Migrate KinematicHardening and its implementations to use typed tensor2 for backstress/flow normals instead of raw arma::vec Voigt vectors. Header changes update signatures for total_backstress, hardening_modulus, update and refresh_state to accept/return tensor2 and clarify strain/stress typing in docs. Implementations (Prager, Armstrong‑Frederick, Chaboche) were adapted to operate on tensor2, use strain/ to_arma_voigt conversions where needed, and fix accumulation/contracting to preserve stress vs. strain Voigt conventions. PlasticityMechanism call sites were updated to convert total_backstress to arma::vec when subtracting from sigma and to pass strain-typed normals into kinematic routines. Tests were updated to expect stress-typed total_backstress and to use strain-typed flow normals. Overall this enforces correct work-conjugacy and Voigt convention handling across kinematic hardening code.
Rename InternalVariable accessors to make raw Voigt/6x6 semantics explicit: vec() -> raw_voigt(), vec_start() -> raw_voigt_start(), mat() -> raw_mat(), mat_start() -> raw_mat_start(). Update header docstrings to explain the "raw" escape-hatch semantics and encourage using set_tensor2()/set_tensor4() for typed sources. Apply corresponding call-site updates across hardening, plasticity, viscoelastic, internal_variable_collection, internal_variable implementation, and unit tests to use the new accessors.
Replace the old requires_rotation API with an explicit is_objective flag to distinguish quantities that are frame-indifferent from rate-like quantities that require custom transport. Scalars are now always objective and their constructor no longer takes a rotate flag. Vector/matrix constructors and InternalVariableCollection registration functions use an objective boolean (instead of rotate). Add typed START accessors (as_tensor2_start, as_tensor4_start) so mechanisms can reconstruct start-of-increment tensorial values safely (used for backward-Euler closed forms). Internals: rename member requires_rotation_ -> is_objective_, update rotate() behaviour to be a no-op for non-objective variables, and add related documentation updates. Refactor InternalVariableCollection to store InternalVariable by-value in a std::deque (registration order == statev layout, stable references) and simplify its API: value-based storage, changed add_* signatures, simplified packing/unpacking/rotation/to_start logic, and removal of several previously-unused helpers. Update callers across mechanisms and tests to the new API and use the new typed start views in hardening closed-form updates.
Move internal-variable ownership into each StrainMechanism and remove the global InternalVariableCollection from ModularUMAT. Mechanisms now own an InternalVariableCollection (ivc_) and provide non-virtual helpers for statev serialization (compute_offsets, pack, unpack, rotate, to_start, variables, statev_size). Registration APIs were simplified to register_variables() with no external ivc argument; many mechanism methods had their InternalVariableCollection parameters removed (compute_constraints, compute_jacobian_contribution, dPhi_dsigma, kappa, inelastic_strain, update, tangent_contribution, compute_work, stiffness_reduction, etc.). ModularUMAT now composes statev as [T_init | mech0 | mech1 | ...], instructs mechanisms to register and compute offsets, and delegates pack/unpack/rotate/to_start per mechanism. Supporting changes: removal of ivc_prefix/key utilities, localized IVC key usage inside mechanisms, small hardening key defaults, and related adjustments across headers and sources to use the new API and ivc_ members. This improves encapsulation, avoids cross-mechanism key collisions, and simplifies hot-loop key handling.
Rename set_solver_params to set_return_mapping_params and clarify comments to indicate these control the local return-mapping (material-point Newton) loop. Remove the older set_max_iterations and set_precision convenience setters. Update tests to call the new API and adjust member comment to "Return-mapping controls".
Pin scikit-build-core to <1 in the GitHub Actions editable install and add --no-build-isolation so the venv's pinned backend builds the editable wheel; include a comment explaining the Windows DLL-vendoring break with 1.x. Also include simcoon/parameter.hpp and replace M_PI with simcoon::pi in a test to avoid MSVC missing M_PI (and add a brief explanatory comment).
Replace the old total_backstress/hardening_modulus API with per-branch compute_backstresses(...) and hardening_modulus(..., X_i) so branch backstresses are built once and passed to the modulus calculation. Add convenience wrappers that sum branch backstresses or compute modulus from ivc for cold paths. Update Prager, Armstrong–Frederick and Chaboche implementations to produce stress-typed branch backstresses and compute the exact tensorial hardening slope (including the <n,n> factor) to avoid Newton overshoots and incorrect tangents. PlasticityMechanism now caches X_branches_, invokes compute_backstresses once per FB iteration and reuses the branches for sigma_eff and hardening modulus. Also remove the conditional gating around applying FB iteration corrections so per-iteration corrections (including negative ds) are always applied; comments explain rationale.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Modular UMAT framework — composable constitutive laws (typed Tensor2/Tensor4 API)
Composable UMAT built from pluggable modules, dispatched as
MODUL(id 200):ElasticityModule— isotropic / cubic / transverse-isotropic / orthotropic; cached typedtensor4stiffness/compliance (L_tensor()/M_tensor()returnconst&, built once per configure).StrainMechanisms —PlasticityMechanism(von Mises/Tresca/Drucker/Hill/DFA/Ani × linear/power-law/Voce/combined-Voce iso × Prager/AF/Chaboche kin),ViscoelasticMechanism(Prony_Nfast generalized Maxwell),DamageMechanism(linear/exponential/power-law/Weibull CDM).ModularUMATorchestrator — convex cutting plane + Fischer–Burmeister multi-constraint solve with 3-phase Jacobian assembly (per-mechanism constraints → cross-mechanism coupling → diagonal); namespaced internal variables (plast0_,visco0_,damage0_, …) with objectivity rotation and statev pack/unpack.simcoon.modular) — frozen dataclasses (ModularMaterial,IsotropicElasticity,Plasticity,Viscoelasticity,Damage, hardening/yield classes) that assembleprops/nstatevforsim.solver("MODUL", ...).Tensor2/Tensor4 integration (this update)
Rebased-by-merge onto
feature/Tensor(Kelvin–Mandel-internaltensor4). The typed layer follows a "typed interface, raw solver core" boundary, chosen by benchmark:dPhi_dsigma()/kappa()returnconst std::vector<tensor2>&, so each operand carries its Voigt convention (dPhi strain-typed; kappa stress-typed for plasticity/viscoelasticity, strain-typed for damage — an inherited convention mix now visible in the type system and documented, not silently altered).dot(a.voigt(), b.voigt())) — bit-identical to the legacy numerics.tensor2fast path for von Mises (Mises/flow_normalon the 3×3, cone-vertex-safe); parametric criteria delegate toEq_stress.plasticity_mechanism.hpp.Bug fixes found while closing test gaps
Bhat(rebuilt locally instead of reusing the cachedK_diag_) → singular Lt, mixed-control solver divergence. Fixed;test_viscoelastic_matches_pronk_referencepins MODUL ≡ reference PRONK to 1.8e-4 relative over a 600-increment relaxation history.StrainMechanism::stiffness_reduction()(damage → 1−D, multiplicative across mechanisms) applied to the elastic prediction; damage integration tests added (elastic uniaxial with unload-freeze, plasticity+damage combined, Python end-to-end).Verification
Tmodular_umat), 6 modular pytests (268 total intest_core), full ctest 41/41.Known modeling items (documented, deliberately not changed here)
DamageMechanism::compute_jacobian_contributionuses a placeholderB(r,r) = 1.0.