Skip to content

Declared reads (3/5): plasticity, hardening, damage and the solvers - #177

Open
petlenz wants to merge 100 commits into
mainfrom
feature/reads-inelastic
Open

petlenz wants to merge 100 commits into
mainfrom
feature/reads-inelastic

Conversation

@petlenz

@petlenz petlenz commented Oct 22, 2026 •

Copy link
Copy Markdown
Member

Part 3 of #168 (epic #173). Stacked on #176.

What

The inelastic materials and the solvers declare their reads:

  • Plasticity:
    • j2_plasticity, j2_rk_plasticity and drucker_prager_plasticity fill stress, tangent, plastic strain and equivalent plastic strain from one callback.
    • They read the strain and, over local edges, the hardening.
    • The two history outputs used to have no callback or writer. They are now written by the return map, so the engine orders them after it.
  • Hardening: linear_isotropic_hardening and exponential_isotropic_hardening read the equivalent plastic strain over a local edge.
  • Damage:
    • isotropic_damage reads its six inputs;
    • scalar_damage_law (power, saturation, Peerlings) fills damage and its derivative from one callback;
    • strain_threshold_yield fills kappa and is_yielding from one callback.
  • autocatalytic_reaction: four callbacks, each with its reads. Residual, jacobian and rate read the history state, so they run after its update, as the creation-order tie-break ran them before.
  • Solvers:
    • backward_euler, rk_integrator: bound with add_callback instead of writing traits().update directly.
    • vector_newton: fills every unknown and every d_<x>_d_strain from its solve. Before, only the first unknown had the callback, and the others were outputs without a writer.

No behaviour change in practice. The return maps, hardening laws and solvers cannot be constant, because they have history outputs or local edges. isotropic_damage and the damage laws now fold when everything they read is constant, which no current deck does. The gain is exact per-callback reads, for #169's edges and #170's check.

The dependencies these solvers hide from the graph (the time step of backward_euler's function, material_ref) are #171.

Tests

No new tests: the plasticity, damage, cure, RK and vector-Newton suites cover these materials unchanged.

Gates (build + ctest incl. ASan variants): GCC 14 1228/1228, g++-13 1228/1228, 0 warnings.

load_driver runs one material point under mixed strain/stress control
(uniaxial stress solves the five lateral strains with the material tangent)
and records the Newton residuals; material_report runs the load cases of a
case.json and writes CSV; reports/make_report.py turns them into figures and
the results section of docs/materials/<material>/README.md.
The correction solves (P:C:P + Isym - P + Iskw) : d_eps = -P:(sigma - sigma*)
with tmech::invf instead of a hand-written 6x6 elimination; Iskw makes the
operator regular on the full 9x9 space. Strain, stress and tangent are
tensors, converted to and from the evaluator buffers with tensor_conversion.
A report document holds the model graph, the material point, recorded
columns wired to any graph property ("material::property[i,j]"), the load
cases, the figures and the checks. Includes merge shared cases and figures;
"$name" strings take parameters. material_report writes resolved.json,
graph.json (materials, parameters, produced and consumed properties) and
CSV per case; make_report.py evaluates the checks against closed forms,
draws the figures and writes the results with the wiring diagram. Every
report is a ctest (checks only).
C++: load_control (stress/strain projectors), load_path (steps as a view),
report_json (parsing into value types), case_runner (material_point and
run_case over a step_sink concept), csv_sink, graph_description; the tool
is now an entry point. load_driver keeps no per-step state: a step_result
carries its start state, so tangent_error works for any earlier step.
Produced properties are listed by name (the registry orders by address).

Python: dataclasses for runs, cases and check results; expressions, checks,
panel renderers (registry by kind) and Markdown sections as separate
functions; unit tests in reports/test_make_report.py (ctest report_script).
The Markdown pages and committed figures are removed. make_report.py
--out <dir> writes <material>.html with embedded figures and an index;
nothing generated is committed. Each report starts with the parameters of
report.json, where each one is used (material, case or check), the derived
values and the material parameters as resolved. resolve_report records
the substitutions in "parameter_usage".
- core/tensor_identities.h: sym/skew identities, volumetric and deviatoric
  projectors, the symmetric 3D basis in Abaqus order; used by
  linear_elasticity, saint_venant_kirchhoff and plasticity_utils
  (make_IIdev / cached_IIdev now delegate; make_isotropic_tangent keeps
  its K I(x)I form).
- core/property_value.h: value of a plain or history property; used by
  material_point_evaluator::resolve_property and the report columns.
- core/graph_description.h: materials with produced (history flag, sorted)
  and consumed properties; the report adds types and parameters.
- umat/load_control.h and umat/material_point_driver.h (was
  report/load_driver.h): mixed strain/stress control of a material point.
make_report.py is the entry point of reports/numsim_report. Expressions
own the document; Raw marks unescaped cells; ReportSummary replaces the
index tuples; axis labels show tensor components as symbols. Tests per
module in reports/tests (ctest report_script discovers them).
parse_point and parse_columns live in report_json.h; step_record and the
step_sink concept in step_record.h, so csv_sink no longer includes the
model builder. load_path::at checks its range, include nesting is capped
at 16 levels, and the finite-difference step of tangent_error is a driver
option (fd_step).
A plain output "iterations", so reports can record the return-map
convergence of the material next to the driver's own Newton.
A check with "plot" draws value and expected per case next to its error
against the tolerance (relative for relative checks); the separate
closed-form figures are gone. Cases take tangent_every (default 1, 0 =
never); the standard cases check every fifth step. The J2 reports record
solver::iterations and show the return map next to the driver's Newton.
…n_property)

Every small-strain consumer took the strain from the fixed property "strain"
of its source. The composites that follow produce one strain per phase
("strain_matrix", "strain_<phase>"), and the phase materials have to read
those. strain_property (default "strain") names the property next to the
existing strain_source / strain_producer_name, the same way tangent_property
and nu_property do; nothing changes for existing documents.

linear_elasticity, linear_stress, j2_plasticity, drucker_prager_plasticity and
the three state functions; tests wire a material to another's stress under
the new name and check the default and a missing property.
"studies" rerun a case with several step counts; the error of a final
quantity against the finest run is fitted with a power of the step size,
and errors below ten times the driver tolerance count as exact. Cases take
"points" (strain states in order) for non-proportional paths; the load
step carries its segment and fraction. The J2 reports show the radial
return exact on a proportional path and first order on tension followed
by shear.
core/elliptic_integrals.h: R_F, R_C, R_D, R_J by Carlson's duplication
algorithm and the complete integrals K, E, so the Eshelby tensors of
ellipsoids and the orientation averaging need no Boost.
core/elastic_constants.h: K, G, E, nu, lambda conversions, the isotropic
stiffness, rotation of rank-2 and rank-4 tensors, rotations from axis and
angle and onto a direction.
umat/deformation_gradient_driver.h drives a material point by F under
mixed F / P control, Newton on the full nine components with dP/dF;
deformation_control names them (21, 31, 32 included). A report whose point
has a deformation_gradient_source records F and P with nine components,
the Cauchy stress as stress_11..23 and J; tensors in cases may have nine
components. Reports for neo_hooke and saint_venant_kirchhoff with exact
checks (S = F^-1 P under uniaxial stress, simple shear, volumetric).
core/eshelby.h: sphere in closed form from the projectors; elliptic and
circular cylinders along any axis or direction (Mura 11.21); ellipsoids
with any semi-axes from Mura's I-integrals through Carlson's R_D, the
degenerate pairs (spheroid, sphere) from the sum rules, an infinite axis
dispatching to the cylinder. Tests against the closed forms, the limits,
uvwmaterial's spheroid formula and the identity S_iikk = (1 + nu)/(1 - nu).
…ses, projectors of the controls, load path accessors
…n tensor, float-safe rotation_to

The elliptic integrals come from Boost.Math where they are first used (the
Eshelby tensors). rotation() builds [axis]_x as -(epsilon . axis); rotation_to()
uses 8 eps of the scalar type instead of 1e-14, which did not compile for float
and would have been c > 1 there. [[nodiscard]] on the pure helpers.
…per shape

eshelby_sphere, eshelby_cylinder, eshelby_spheroid and eshelby_ellipsoid consume
nu_source / nu_property and publish the output "tensor"; the engine orders and
re-evaluates them with nu. Two scalar producers feed them: poisson_ratio_from_moduli
(K_source, G_source -> value) and isotropic_moduli (tangent_source -> K, G,
poisson_ratio through isotropize, rejecting an anisotropic tangent). Composites
hold a node through eshelby_tensor_interface (tensor() at the wired nu, tensor(nu)
for the self-consistent iterate). core/eshelby.h keeps the closed forms: the
integrals come from Boost.Math (ellint_rd), the tensors are assembled on the
eigenprojectors, equal axes are detected at sqrt(eps).
core/isotropize.h: K = Pvol::C / Dim, G = Pdev::C / (2 (Dim(Dim+1)/2 - 1)),
and the isotropic part built from them. core/lebedev.h with generated
tables (tools/extract_lebedev.py from uvwmaterial's rule, precisions 3 to
23, up to 194 points): spherical means and integrals of tensor-valued
functions of the direction. eshelby::depolarisation for ellipsoids via
Carlson's R_J, the elliptic cylinder in closed form, and the cylinder along
a direction.
Phases are materials providing a tangent, with volume fractions; outputs
the effective tangent, the stress and the strain concentration tensors
(sum c_i A_i = Isym). Report with the closed forms of the three kernels
and their ordering. Report generator: derived values resolve in
dependency order, independent of the key order.
…ent producer over phases

The composite consumes the phases' produced tangents, fractions and Eshelby
tensors and produces the effective tangent and the concentrations; the stress
is a linear_stress on the tangent. It reaches into no other node: no
material_ref, no eshelby(nu), no isotropisation of the matrix.

Self-consistent is a graph cycle closed by one local edge: isotropic_moduli
on the composite tangent produces nu of the effective medium, the Eshelby
nodes read it with nu_edge = local (from nu_initial on the first pass), and
one update() is one relaxed Picard step with "residual" as its error. The
matrix is a phase with its own Eshelby node, so its shape is no longer a
hard-coded sphere. material_point_evaluator closes the fixed point per host
call (residual_property, residual_tolerance, max_passes; the report point
names them). isotropic_moduli with a non-positive tolerance projects without
the isotropy check, which the effective medium needs.
…ution

Spheroidal fibres in an isotropic matrix with the second-order orientation
tensor: Mori-Tanaka with the orientation-averaged dilute concentration, or
the averaged unidirectional stiffness. Eshelby tensors per direction are
cached for the matrix Poisson ratio. Report with aligned, injection-moulded
and random states against the unidirectional composite and isotropy.
… phase, passes to the fixed point

A phase whose tangent depends on its own strain closes a cycle through
mean_field_strain: the concentrations need the algorithmic tangent, the
tangent needs the local strain. The phase node reads that tangent with
tangent_edge = local (no ordering, the value of the previous pass; the
initial_tangent_source on the first pass), so one update() is one pass and
mean_field_stress::residual says when the passes agree. The evaluator
closes the fixed point per host call through residual_property.

Tests with an isotropic-damage matrix and stiff spheres: below the threshold
the chain equals the linear composite; after damage the local strains
partition the composite strain, the state function saw the matrix strain, the
stress is the phase average, the tangent is the concentration-weighted sum and
matches central differences on a vanishing step; the residual decreases over
the passes of a large step; stress control through the driver converges; the
Reuss kernel ends with equal stresses in both phases. The report
mean_field_damage runs the standard cases with the point closing the fixed
point; the card converter maps a rule-of-mixture card with a nonlinear phase
to the pair for every kernel.
core/log_strain.h: ln V and d ln V / dF from an eigendecomposition of
F F^T (Eigen) with the divided-difference derivative, continuous through
repeated eigenvalues; tmech's log of V gives wrong values at eigenvalue
splits of about 1e-7 (petlenz/tmech#37). hencky_elasticity: tau = C :
(ln V - eps_eig), P = tau F^-T, consistent dP/dF, optional eigenstrain.
Report with the exact uniaxial and volumetric closed forms; reports take
a per-case fd_step.
docs/uvwmaterial-mapping.md: every class of uvwmaterial with its
numsim-materials counterpart, parameters and status; the card formats.
tools/convert_uvwmaterial_card.py turns the nested cards of the
supported types into json_model documents, naming the issue of the port
for the others; tested on the example cards, and the converted
Mori-Tanaka card loads and gives the analytic moduli.
…se with a hardening J2 matrix

mean_field_damage records the intact Mori-Tanaka next to the local pair
(composite and phase stresses, damage, local strains under strain and stress
control; on first loading the damaged composite is softer, after a reversal
the incremental localisation leaves a residual stress, noted). The new case
mean_field_reuss_plastic runs the Reuss kernel through the passes with a J2
matrix: the phases carry the same stress to first order in the step.
…he phase materials

The composite is split into two nodes. mean_field_strain consumes the composite
strain increment and the phases' tangents, fractions and Eshelby tensors and
produces one strain per phase (history) and the concentration tensors, for
the voigt, reuss, hill, dilute and mori_tanaka kernels. The phase materials
consume those strains through strain_source / strain_property like any other
strain. mean_field_stress averages the stresses and tangents they produce with
the concentrations and reports a residual, the difference between the tangent
a phase node carried and its material's current one.

With linear phases every edge is global, one pass is exact and the result
equals mean_field_tangent with a linear_stress to round-off, the phase strains
and stresses available as properties. docs/mean_field.md describes the
composite layer.
tmech main (46d8fdd, after tmech#38) is exact at repeated eigenvalues, so the
material uses ln b from tmech::log(b) and its derivative, is Dim generic, and
core/log_strain.h with its Eigen dependency goes away. The per-case fd_step of
the reports is optional and no longer overrides the driver default; the 1e-6
workaround in the finite-strain standard is dropped (all three finite-strain
reports pass at 1e-7). New tests: tangent at repeated and nearly repeated
stretches (the tmech#37 symptom), objectivity P(RU) = R P(U) and P(R) = 0.
Mapping document and card converter list the Hencky family.
Main has the reviewed split of this branch (#199-#203), so main's version is
taken everywhere; the branch no longer differs from main.
…ls-and-constants

# Conflicts:
#	tests/CMakeLists.txt
@petlenz
petlenz changed the base branch from feature/reads-elasticity to main October 24, 2026 21:17
…eshelby-tensors

# Conflicts:
#	tests/CMakeLists.txt
…v-depolarisation

# Conflicts:
#	tests/CMakeLists.txt
…/rule-of-mixture

# Conflicts:
#	tests/CMakeLists.txt
# Conflicts:
#	include/numsim-materials/report/report_json.h
#	tests/CMakeLists.txt
residual_source, residual_tolerance and max_passes are driver parameters
(material_point's defaults) and keys of the report's "point"; the mean-field
tangent report uses them for the self-consistent fixed point.
…averaging

# Conflicts:
#	tests/CMakeLists.txt
…-composite

# Conflicts:
#	tests/CMakeLists.txt
…-elasticity

# Conflicts:
#	include/numsim-materials/umat/deformation_gradient_driver.h
#	include/numsim-materials/umat/material_point_driver.h
#	tests/CMakeLists.txt
…genstrain

# Conflicts:
#	tests/CMakeLists.txt
# Conflicts:
#	include/numsim-materials/report/report_json.h
#	tests/CMakeLists.txt
#	tests/test_report_case.cpp
The report's "point" names them for materials that publish e.g.
damaged_stress / damaged_tangent (isotropic_damage report).
# Conflicts:
#	reports/CMakeLists.txt
#	tests/CMakeLists.txt
# Conflicts:
#	include/numsim-materials/report/case_runner.h
#	include/numsim-materials/report/report_json.h
#	tests/CMakeLists.txt
…eld-composite-global

# Conflicts:
#	docs/materials/mean_field_composite_global/report.json
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant