Repository navigation
Isotropisation, Lebedev quadrature and depolarisation tensors - #129
Merged
Merged
Conversation
This was referenced Oct 16, 2026
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 16, 2026 23:43
b33031f to
e0f31b7
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 16, 2026 23:43
782e7f9 to
94b5a8a
Compare
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. isotropize / Lebedev / depolarisation — needs changes (small)
|
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 17, 2026 13:31
e0f31b7 to
f949af6
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 17, 2026 13:37
94b5a8a to
a9613ff
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 17, 2026 13:39
f949af6 to
0747e40
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 17, 2026 13:45
a9613ff to
0050150
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 17, 2026 15:29
0747e40 to
04688cb
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 17, 2026 15:30
0050150 to
191fe75
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 17, 2026 21:02
04688cb to
7fad660
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 17, 2026 21:07
191fe75 to
dea9c37
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
2 times, most recently
from
October 17, 2026 22:20
12e4373 to
1d5e486
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 17, 2026 22:25
dea9c37 to
1141085
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 18, 2026 11:27
1d5e486 to
6b42365
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 18, 2026 11:34
1141085 to
daf4cd1
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 18, 2026 13:38
6b42365 to
1dfcb1c
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 18, 2026 13:38
daf4cd1 to
361a892
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 18, 2026 20:03
1dfcb1c to
270489d
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 18, 2026 20:03
361a892 to
6e9cbdc
Compare
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 18, 2026 20:22
270489d to
3c9839a
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 18, 2026 20:22
6e9cbdc to
1d143ff
Compare
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.
petlenz
force-pushed
the
feature/eshelby-tensors
branch
from
October 19, 2026 04:06
3c9839a to
61865b0
Compare
petlenz
force-pushed
the
feature/isotropize-lebedev-depolarisation
branch
from
October 19, 2026 04:06
1d143ff to
207536a
Compare
…ls-and-constants # Conflicts: # tests/CMakeLists.txt
…eshelby-tensors # Conflicts: # tests/CMakeLists.txt
…v-depolarisation # Conflicts: # tests/CMakeLists.txt
petlenz
marked this pull request as ready for review
October 27, 2026 17:49
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
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.
Closes #100, closes #101, closes #102. Part of the uvwmaterial port (#94); stacked on #128.
core/isotropize.h:isotropize<T, Dim>(C)→{K, G}of the closest isotropic tensor,K = Pvol::C / Dim,G = Pdev::C / (2 (Dim(Dim+1)/2 − 1))(the formulas oftensor_isotropization_general, written for any dimension), andisotropic_part(C). Tests: an isotropic C is its own projection; a deviatoric rank-one term shifts only G, by the predicted amount; the moduli are rotation invariant for an anisotropic C.core/lebedev.hwithcore/lebedev_tables.hgenerated bytools/extract_lebedev.pyfrom uvwmaterial'ssphere_lebedev_rule(Lebedev–Laikov): precisions 3, 5, 7, 9, 11, 13, 15, 23 (6 to 194 points) as unit vectors and weights summing to one, instead of the 13.9k-line table.rule(p),average(p, f)andintegrate(p, f)for scalar or tensor-valuedf(n). Tests: weights sum to one and nodes are unit; the second moment is I/3 and the fourth (I⊗I + 2 Isym)/15 for every precision; the degree of exactness (the 6-point rule misses n₁⁴, the 14-point rule gets it).eshelby::depolarisation(a1, a2, a3): N_ii = (a₁a₂a₃/3) R_J(a₁², a₂², a₃², a_i²); an infinite axis gives the elliptic cylinder N_i = a_j/(a_i + a_j);depolarisation_cylinder(direction). The original's elliptic-cylinder formula (a spheroid expression) is not carried over. Tests: sphere I/3, needle and disc limits, trace one, the elliptic cylinder, prolate and oblate spheroids against the closed forms (Osborn).818/818 tests pass (GCC 14, reports on).
JSON input
New material types:
depolarisation_sphere,depolarisation_cylinder,depolarisation_spheroid,depolarisation_ellipsoid. Ajson_modeldocument that loads (theFromJsontest of the PR):{"materials": [ {"type": "depolarisation_sphere", "name": "sphere"}, {"type": "depolarisation_cylinder", "name": "cylinder", "direction": [0, 1, 0]}, {"type": "depolarisation_spheroid", "name": "spheroid", "aspect_ratio": 5, "direction": [1, 0, 0]}, {"type": "depolarisation_ellipsoid", "name": "ellipsoid", "semi_axes": [1, 2, 3]} ]}Producers without inputs (the conductivity Eshelby tensor of an isotropic matrix is geometry only); output
tensor(rank 2, trace one), and nothing else — no interface, composites consume the property. Aphase(#130) republishes it through"depolarisation_source": "sphere"for the conductivity composites (#110).