Skip to content

inclusion and mean_field_composite: dilute, Mori-Tanaka and self-consistent - #131

Merged
petlenz merged 30 commits into
mainfrom
feature/mean-field-composite
Oct 27, 2026
Merged

petlenz merged 30 commits into
mainfrom
feature/mean-field-composite

Conversation

@petlenz

@petlenz petlenz commented Oct 16, 2026 •

Copy link
Copy Markdown
Member

Closes #107, closes #143. Part of the uvwmaterial port (#94); stacked on #130.

mean_field_tangent: dilute, Mori–Tanaka and self-consistent effective stiffness as a tangent producer over phase nodes (#130). matrix and phases are phases (produced tangent, fraction and Eshelby tensor); kernel dilute | mori_tanaka | self_consistent. Outputs tangent, concentration_<phase> for every phase (the matrix included) and residual; the stress is a linear_stress on the tangent. Dilute: Tᵢ = (I + Sᵢ : Cₘ⁻¹ : (Cᵢ − Cₘ))⁻¹; Mori–Tanaka: Aᵢ = Tᵢ : (cₘ I + Σ cⱼ Tⱼ)⁻¹. The Eshelby tensors come in as properties, so their Poisson's ratio is whatever node produces it (isotropic_moduli of the matrix tangent, or a constant) — the composite reaches into no other node and isotropises nothing.

Self-consistent as a graph cycle. Every phase, the matrix included (its own Eshelby node — closes #143), is embedded in the effective medium, whose ν is isotropic_moduli on this node's tangent (tolerance: 0, the projection without the isotropy check). The Eshelby nodes consume that ν over a local edge ("nu_edge": "local", starting from nu_initial): EdgeKind::Local imposes no ordering, so one update() is one relaxed Picard step (relaxation, default 0.5) on the tangent of the previous pass, and residual = ‖C_new − C‖/‖C_new‖. The fixed point is run by whoever calls update():

  • material_point_evaluator gets residual_property, residual_tolerance, max_passes and repeats the graph update per host call until the residual is below tolerance (report documents: "point": {..., "residual": "scs::residual", "tolerance": 1e-10, "max_passes": 200}).
  • A solver's Newton loop re-evaluates the material anyway; the driver and the UMAT close the point through the evaluator.

The old version iterated inside the node with material_ref calls eshelby(nu) into the inclusions and a hard-coded sphere for the matrix; both are gone. isotropic_moduli with a non-positive tolerance projects without the check (needed for the effective medium).

Tests (test_mean_field_tangent, 13): a phase republishes sphere, cylinder and spheroid tensors from its node; a phase without a shape is rejected at finalize(); Mori–Tanaka spheres match the analytic K and G (= Hashin–Shtrikman lower bound) at c = 0.05, 0.3, 0.6 to 1e-10 with zero residual; dilute matches its closed-form K and the Mori–Tanaka limit as c → 0; two phases of one material equal one; voids soften with the exact K; Σ c A = I and the linear_stress gives C : ε for every kernel; the self-consistent fixed point is symmetric under swapping matrix and inclusion (1e-8) and stiffer than Mori–Tanaka, converges within the pass loop and stays put on another pass; the evaluator closes the fixed point in one call; aligned fibres are transversely isotropic with E_L near Voigt; overfull composites and unknown kernels throw; JSON.

Report docs/materials/mean_field_tangent/report.json: Mori–Tanaka, dilute, self-consistent (its own Eshelby node at ν_eff over the local edge, converged per step through the point's residual) and unidirectional on one point, the closed forms as checks.

862/862 tests pass (GCC 14, reports on).

JSON input

New material type: mean_field_tangent. A json_model document that loads (the FromJson test of the PR):

{"materials": [
  {"type": "external_strain_source", "name": "strain"},
  {"type": "constant_scalar", "name": "Km", "value": 3.0},  {"type": "constant_scalar", "name": "Gm", "value": 1.0},
  {"type": "constant_scalar", "name": "Ki", "value": 60.0}, {"type": "constant_scalar", "name": "Gi", "value": 30.0},
  {"type": "isotropic_tangent", "name": "Cm", "K_source": "Km", "G_source": "Gm"},
  {"type": "isotropic_tangent", "name": "Ci", "K_source": "Ki", "G_source": "Gi"},
  {"type": "isotropic_moduli", "name": "moduli", "tangent_source": "Cm"},
  {"type": "eshelby_sphere", "name": "sphere", "nu_source": "moduli", "nu_property": "poisson_ratio"},
  {"type": "constant_scalar", "name": "c", "value": 0.3},
  {"type": "constant_scalar", "name": "cm", "value": 0.7},
  {"type": "phase", "name": "spheres", "tangent_source": "Ci", "fraction_source": "c", "eshelby_source": "sphere"},
  {"type": "phase", "name": "matrix", "tangent_source": "Cm", "fraction_source": "cm"},
  {"type": "mean_field_tangent", "name": "mf", "matrix": "matrix", "phases": ["spheres"], "kernel": "mori_tanaka"},
  {"type": "linear_stress", "name": "composite", "tangent_source": "mf", "strain_source": "strain"}
]}

Self-consistent: the matrix phase carries an Eshelby node too, both nodes read ν of the effective medium over a local edge, and the point closes the fixed point:

{"materials": [
  {"type": "isotropic_moduli", "name": "moduli_eff", "tangent_source": "mf", "tolerance": 0},
  {"type": "eshelby_sphere", "name": "sphere", "nu_source": "moduli_eff", "nu_property": "poisson_ratio", "nu_edge": "local", "nu_initial": 0.35},
  {"type": "phase", "name": "spheres", "tangent_source": "Ci", "fraction_source": "c", "eshelby_source": "sphere"},
  {"type": "phase", "name": "matrix", "tangent_source": "Cm", "fraction_source": "cm", "eshelby_source": "sphere"},
  {"type": "mean_field_tangent", "name": "mf", "matrix": "matrix", "phases": ["spheres"], "kernel": "self_consistent", "relaxation": 0.5},
  {"type": "linear_stress", "name": "composite", "tangent_source": "mf", "strain_source": "strain"}
 ],
 "point": {"strain_source": "strain", "stress_source": "composite", "tangent_source": "mf", "residual": "mf::residual", "tolerance": 1e-10, "max_passes": 200}}

Outputs tangent, concentration_<phase> (matrix included), residual. Several inclusions: more phase entries (each with its own Eshelby node) in phases.

Update after main (#199–#205)

The report's "point" now builds the driver's parameters, so this PR carries #210: the drivers take residual_source, residual_tolerance and max_passes (material_point's defaults) and pass them to the material point, and "point" accepts the same keys. docs/materials/mean_field_tangent/report.json now uses "residual_source": "scs::residual", "residual_tolerance": 1e-10 instead of "residual" / "tolerance" (which is the driver's Newton tolerance on main). New tests: MaterialPointDriver.PassesTheResidualSettingsToThePoint, DeformationGradientDriver.PassesTheResidualSettingsToThePoint, and the residual keys and max_passes >= 1 in test_report_case.

Closes #210.

@petlenz

petlenz commented Oct 17, 2026

Copy link
Copy Markdown
Member Author

Review findings (design, tmech expressions instead of component loops, modern C++, comments, tests, correctness). Severity: must = wrong result, hang, build break or unreachable feature; should = design, duplication, accuracy; nit = style.

inclusion + mean_field_composite — needs changes

  • [must] SC fixed point with isotropised Eshelby tensors and mixed shapes → tangent without major symmetry: ‖C−Cᵀ‖/‖C‖ = 4.6e-2 (spheres 0.2 + spheroids ar 10, probe); MT/dilute ≤ 3e-16. Symmetrise each iterate and the result, or restrict SC to isotropic media and throw; add a major-symmetry test for mixed shapes.
  • [should] C = (C + C_new)/2 unconditional: with normalised concentrations plain iteration converges in 17–27 iterations, relaxed in 46–67 (probe). Drop or expose relaxation (default 1). The divergence cited in the PR body was uvwmaterial's unnormalised kernel.
  • [should] dilute sets A_matrix = Is (Σ c A ≠ I, test special-cases it) → (Is − Σ cᵢTᵢ)/cm as short_fibre does.
  • [should] :111-116 duplicates MT start verbatim; poisson_of duplicated in short_fibre → isotropize.h: require_isotropic(C, tol); std::clamp(poisson, −0.99, 0.49) hides an iterate leaving the admissible set; Eshelby recomputed every update (memoise in inclusion::eshelby(ν)).
  • [should] inclusion.h:53-54 direction copy loop → tmech::adaptor<value_type const,3,1,tmech::full<3>>(d.data()); volume_fraction required even with fraction_source (tests write a dummy 0.0).
  • [should] static_assert(Dim == 3) inside types registered unconditionally → register_default_materials<material_policy_2d>() breaks; register under if constexpr (Dim == 3).
  • [nit] convergence tested after relaxation; no non-convergence test; ellipsoid never built through inclusion; SC only checked by ordering/swap (Hill SC closed form for spheres would pin 1e-10).

@petlenz

petlenz commented Oct 17, 2026

Copy link
Copy Markdown
Member Author

Rebased on the Eshelby-node refactor of #128: inclusion now names its eshelby node instead of carrying shape parameters (the uvwmaterial model, one tensor object per inclusion). Of the review findings this commit also covers:

  • should, poisson_of duplicated: gone, require_isotropic (core/isotropize.h) for the matrix, isotropize + poisson_ratio for the self-consistent iterate.
  • should, Eshelby recomputed every update: the node caches its tensor by ν and eshelby(nu) returns the cached value when ν is the matrix's.
  • should, direction copy loop: moved into the node (tmech::adaptor<…, full<3>>).
  • should, static_assert(Dim == 3) in unconditionally registered types: inclusion and mean_field_composite register under if constexpr (Dim == 3).

Still open from the review (not part of the refactor): symmetrising the self-consistent tangent for mixed shapes, the ½ relaxation, dilute's A_matrix, the Mori–Tanaka start duplicated as the SC start, volume_fraction required beside fraction_source (#137).

@petlenz

petlenz commented Oct 17, 2026

Copy link
Copy Markdown
Member Author

Producer/consumer pass: inclusion is a phase of any composite — material, fraction_source (the only way to give a fraction; the volume_fraction parameter is gone, so the parameter-or-source duality the review flagged cannot occur) and an optional eshelby node. The Eshelby nodes themselves consume nu_source/nu_property (#128): fixtures and the report wire isotropic_moduli of the matrix tangent to them, so a matrix whose tangent evolves re-evaluates its Eshelby tensors through the engine. mean_field_composite rejects a phase without a shape at update.

@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from 7d5a01f to 1a2e21a Compare October 17, 2026 21:07
@petlenz
petlenz force-pushed the feature/mean-field-composite branch from 9c83c8c to 2bcebfa Compare October 17, 2026 21:13
@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from 1a2e21a to 13a8bd0 Compare October 17, 2026 22:32
@petlenz
petlenz force-pushed the feature/mean-field-composite branch from 2bcebfa to c88bed0 Compare October 17, 2026 22:38
@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from 13a8bd0 to f64e188 Compare October 18, 2026 11:42
@petlenz
petlenz force-pushed the feature/mean-field-composite branch 2 times, most recently from 5434109 to 7bd385b Compare October 18, 2026 12:02
@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from f64e188 to 6528f13 Compare October 18, 2026 13:38
@petlenz
petlenz force-pushed the feature/mean-field-composite branch from 7bd385b to c23516a Compare October 18, 2026 13:38
@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from 6528f13 to 5d53585 Compare October 18, 2026 20:03
@petlenz
petlenz force-pushed the feature/mean-field-composite branch from c23516a to 4e04cb7 Compare October 18, 2026 20:03
@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from 5d53585 to 57ada14 Compare October 18, 2026 20:22
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.
@petlenz
petlenz force-pushed the feature/rule-of-mixture branch from 57ada14 to 6520437 Compare October 19, 2026 04:06
@petlenz
petlenz force-pushed the feature/mean-field-composite branch from c43505c to 655efa4 Compare October 19, 2026 04:06
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/rule-of-mixture to main October 24, 2026 21:16
…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.
@petlenz
petlenz marked this pull request as ready for review October 27, 2026 17:49
@petlenz
petlenz merged commit 2dfc16d into main Oct 27, 2026
2 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

1 participant