Skip to content

Eshelby tensors as materials: sphere, cylinder, spheroid and ellipsoid nodes - #128

Merged
petlenz merged 23 commits into
mainfrom
feature/eshelby-tensors
Oct 27, 2026
Merged

petlenz merged 23 commits into
mainfrom
feature/eshelby-tensors

Conversation

@petlenz

@petlenz petlenz commented Oct 16, 2026 •

Copy link
Copy Markdown
Member

Closes #96, closes #97, closes #98, closes #99, closes #103. Part of the uvwmaterial port (#94); stacked on #127 (which sits on #145).

Each Eshelby tensor is a material, and a pure producer. As in uvwmaterial (eshelby_tensor_solid_sphere, _cylinder_x1, _cylindrical_ellipsoid, _general_geometry) every shape is its own node of the graph:

  • eshelby_sphere, eshelby_cylinder (direction), eshelby_spheroid (aspect_ratio = length over diameter, direction), eshelby_ellipsoid (semi_axes along the global axes).
  • Every node consumes Poisson's ratio (nu_source, nu_property, default value) and publishes the tensor as the output tensor. Nothing else: no interface, no tensor(nu) for callers off the graph. A composite consumes tensor like any other property; the engine re-evaluates the node when ν changes, so a curing matrix is followed and a constant one costs one evaluation. The common plumbing lives once in the CRTP detail::eshelby_tensor_base; a shape is ~20 lines (the callback is bound from the most-derived constructor, UBSan-clean).
  • Registered for Dim = 3 only (if constexpr), so a 2D registration still compiles.

Producers of ν: poisson_ratio_from_moduli takes K and G as parameters and computes ν once at construction (output value, no callback, like constant_scalar) — the constant case. isotropic_moduli consumes a tangent (tangent_source, tangent_property) and produces K, G, poisson_ratio per update, rejecting a tangent off its isotropic projection by more than tolerance — the case of a ν that follows a tangent.

core/eshelby.h keeps the closed forms (Mura ch. 11) as free functions behind the nodes: sphere(nu) from the projectors, elliptic_cylinder(nu, a, b, axis), cylinder(nu, axis | direction), ellipsoid(nu, a1, a2, a3) with I_i = (4π/3) a₁a₂a₃ R_D from Boost.Math (ellint_rd; an installed Boost ≥ 1.74 or the standalone boostorg/math repository fetched at boost-1.83.0 — a header-only dependency; like Eigen it is exported only when installed). The tensors are assembled on the eigenprojectors eᵢ⊗eᵢ (detail::orthotropic(A, s) with otimes/otimesu/otimesl) instead of per component, and equal axes are detected at √ε so the difference rule (I_j − I_i)/(a_i² − a_j²) never runs into its cancellation window (continuity into the spheroid tested at δ = 1e-9, 1e-11, 1e-13).

Tests: test_eshelby (8, the closed forms); test_eshelby_materials (8, the nodes: each shape against its closed form at the consumed ν, a changing ν re-evaluates the node, ν from an isotropic tangent through isotropic_moduli, invalid geometry rejected at construction, all four from JSON); test_scalar_producers (4: ν fixed at construction, non-positive moduli rejected, the K-G-tangent-moduli round trip, JSON); test_isotropize (5).

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

JSON input

New material types: poisson_ratio_from_moduli, isotropic_moduli, eshelby_sphere, eshelby_cylinder, eshelby_spheroid, eshelby_ellipsoid. A json_model document that loads (the FromJson test of the PR):

{"materials": [
  {"type": "poisson_ratio_from_moduli", "name": "nu", "K": 100, "G": 40},
  {"type": "eshelby_sphere", "name": "sphere", "nu_source": "nu"},
  {"type": "eshelby_cylinder", "name": "cylinder", "nu_source": "nu", "direction": [0, 0, 1]},
  {"type": "eshelby_spheroid", "name": "spheroid", "nu_source": "nu", "aspect_ratio": 10, "direction": [1, 0, 0]},
  {"type": "eshelby_ellipsoid", "name": "ellipsoid", "nu_source": "nu", "semi_axes": [1, 2, 3]}
]}

When the matrix tangent is itself a product (cure, temperature, a composite), isotropic_moduli turns it into K, G and ν; Eshelby nodes then consume "nu_source": "moduli", "nu_property": "poisson_ratio":

{"materials": [
  {"type": "constant_scalar", "name": "K", "value": 3333.33},
  {"type": "constant_scalar", "name": "G", "value": 1153.85},
  {"type": "isotropic_tangent", "name": "stiffness", "K_source": "K", "G_source": "G"},
  {"type": "isotropic_moduli", "name": "moduli", "tangent_source": "stiffness"},
  {"type": "eshelby_sphere", "name": "sphere", "nu_source": "moduli", "nu_property": "poisson_ratio"}
]}

Outputs: value (ν), K/G/poisson_ratio, tensor (rank 4).

@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.

Eshelby — needs changes

  • [should] eshelby.h:118 equal threshold 1e-12 leaves a cancellation window: ellipsoid(ν,1,1+δ,2) vs spheroid 3.9e-5 at δ=1e-12, 8.6e-7 at 1e-10. Use sqrt(eps); test δ=1e-11.
  • [should] :29-45, :66-78 per-index assembly, shear block repeated four times → S = Σ A_ij otimes(E_i,E_j) + Σ_{i≠j} s_ij (otimesu(E_i,E_j)+otimesl(E_i,E_j)) with E_i = otimes(e_i,e_i); one detail::orthotropic(A, s) serves from_integrals and elliptic_cylinder.
  • [should] semi-axis validation / infinite-axis parsing duplicated (:98-108, :139-148) → detail::semi_axes; std::size_t axis + run-time range check + ambiguous overload pair → enum class axis.
  • [should] header must state: ν is the matrix Poisson ratio; S has minor symmetries only.
  • [nit] cyclic convention of elliptic_cylinder undocumented; <limits> unused; test helper trace computes S_iikk; spheroid reference is uvwmaterial's own transcription (cite Mura 11.28/11.29); no ν<0, no float.

@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from e0f31b7 to f949af6 Compare October 17, 2026 13:31
@petlenz petlenz changed the title Eshelby tensors: sphere, cylinders, spheroids and general ellipsoids Eshelby tensors as materials: sphere, cylinder, spheroid and ellipsoid nodes Oct 17, 2026
@petlenz

petlenz commented Oct 17, 2026

Copy link
Copy Markdown
Member Author

Refactored as requested: each Eshelby tensor is its own material (eshelby_sphere, eshelby_cylinder, eshelby_spheroid, eshelby_ellipsoid), wired to the matrix by name and publishing the output tensor; the PR body describes the design. The consumers (inclusion, mean_field_composite in #131, short_fibre_composite in #133, the curing test in #137, the converter in #139) follow as the stack is rebased.

Review findings addressed in the same commit:

  • should, per-index assembly: detail::orthotropic(A, s) on the eigenprojectors, shared by the ellipsoid and the elliptic cylinder.
  • should, near-equal axes: threshold √ε; continuity tested at δ = 1e-9, 1e-11, 1e-13 to 1e-8.
  • should, semi-axis parsing duplicated: detail::infinite_axis.
  • should, header doc: ν is the matrix Poisson's ratio, minor symmetries only.
  • Carlson integrals: Boost.Math ellint_rd (dependency added here, where it is first used).

Kept: std::size_t axis for the free cylinder functions (now internal to the nodes); the depolarisation tensors stay in #129.

@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from f949af6 to 0747e40 Compare October 17, 2026 13:39
@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from 0747e40 to 04688cb Compare October 17, 2026 15:29
@petlenz
petlenz force-pushed the feature/elliptic-integrals-and-constants branch from 17946ba to 4313828 Compare October 17, 2026 20:57
@petlenz
petlenz force-pushed the feature/eshelby-tensors branch 4 times, most recently from 1d5e486 to 6b42365 Compare October 18, 2026 11:27
@petlenz petlenz closed this Oct 18, 2026
@petlenz petlenz reopened this Oct 18, 2026
@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from 6b42365 to 1dfcb1c Compare October 18, 2026 13:38
@petlenz
petlenz force-pushed the feature/elliptic-integrals-and-constants branch from 4313828 to bb455b7 Compare October 18, 2026 13:41
@petlenz petlenz closed this Oct 18, 2026
@petlenz petlenz reopened this Oct 18, 2026
@petlenz
petlenz force-pushed the feature/elliptic-integrals-and-constants branch from bb455b7 to 56e0364 Compare October 18, 2026 20:03
@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from 1dfcb1c to 270489d Compare October 18, 2026 20:03
@petlenz
petlenz force-pushed the feature/elliptic-integrals-and-constants branch from 56e0364 to 25c92c6 Compare October 18, 2026 20:22
@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from 270489d to 3c9839a Compare October 18, 2026 20:22
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).
@petlenz
petlenz force-pushed the feature/elliptic-integrals-and-constants branch from 25c92c6 to 405ca6c Compare October 19, 2026 04:06
@petlenz
petlenz force-pushed the feature/eshelby-tensors branch from 3c9839a to 61865b0 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/elliptic-integrals-and-constants to main October 24, 2026 21:16
…eshelby-tensors

# Conflicts:
#	tests/CMakeLists.txt
@petlenz
petlenz marked this pull request as ready for review October 27, 2026 17:49
@petlenz
petlenz merged commit a0c8fa1 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