API — Localization & contribution

strain_strain_loc(incl, C₁, C₀; kw...)  -> Tens{4,3}

Dilute strain-strain localization tensor A_εε (Eshelby).

stress_strain_loc(incl, C₁, C₀; kw...)  -> Tens{4,3}

strain_stress_loc(incl, C₁, C₀; kw...)  -> Tens{4,3}

stress_stress_loc(incl, C₁, C₀; kw...)  -> Tens{4,3}

gradient_gradient_loc(incl, K₁, K₀; kw...) -> Tens{2,3}

flux_gradient_loc(incl, K₁, K₀; kw...)    -> Tens{2,3}

gradient_flux_loc(incl, K₁, K₀; kw...)    -> Tens{2,3}

flux_flux_loc(incl, K₁, K₀; kw...)        -> Tens{2,3}

stiffness_contribution(incl, C₁, C₀; kw...) -> Tens{4,3}
stiffness_contribution(crack, C₀; kw...)   -> Tens{4,3}

Size-independent stiffness contribution tensor N of an inclusion in a matrix C₀. For a dilute family of volume fraction f: ΔC_eff = f · N (see delta_stiffness).

compliance_contribution(incl, P₁, P₀; kw...) -> Tens
compliance_contribution(incl, P₀; kw...)     -> Tens

Size-independent compliance contribution tensor H of an inclusion in a matrix P₀. For a dilute family of volume fraction f: ΔS_eff = f · H (see delta_compliance).

The two-argument form is the flat-object flavour used by cracks and, more generally, by any inclusion registered in an RVE with a CrackDensity amount: the inclusion carries no property of its own, and the amount is applied afterwards by the three-argument delta_compliance.

Methods for the built-in cracks live in Cracks/compliance.jl; the generic three-argument method for solid inclusions lives in src/contribution.jl.

conductivity_contribution(incl, K₁, K₀; kw...) -> Tens{2,3}
conductivity_contribution(crack, K₀; kw...)     -> Tens{2,3}

Size-independent conductivity contribution tensor for the 2nd-order transport problem. Analog of stiffness_contribution.

resistivity_contribution(incl, K₁, K₀; kw...) -> Tens{2,3}

Size-independent resistivity contribution tensor of an inclusion (2nd-order analog of compliance_contribution for solid ellipsoids).

delta_stiffness(N, f) -> Tens{4,3}

Dilute effective-stiffness correction ΔC = f · N from the size- independent contribution tensor N and the volume fraction f.

delta_compliance(H, f)        -> Tens
delta_compliance(incl, H, ε)  -> Tens

Dilute effective-compliance correction from the size-independent contribution tensor H.

The two-argument form is the volume-fraction one, ΔS = f · H. The three-argument form is the amount × contribution seam of flat objects: it carries the geometric prefactor relating a density-like amount to the effective correction (4π/3 for an elliptical crack of Budiansky density ε³ᵈ = N a b², π for a ribbon crack of ε²ᵈ = N b²). Every inclusion meant to be registered with a CrackDensity amount must provide the three-argument methods of delta_compliance, delta_stiffness, delta_conductivity and delta_resistivity.

delta_conductivity(N_K, f) -> Tens{2,3}

Dilute effective-conductivity correction ΔK = f · N_K.

delta_resistivity(R, f)        -> Tens{2,3}
delta_resistivity(incl, R, ε)  -> Tens{2,3}

Dilute effective-resistivity correction — 2nd-order analog of delta_compliance, with the same two calling conventions.

is_homogeneous_inclusion(incl) -> Bool

Whether the inclusion carries a single uniform property, so that the mean-field identities of a homogeneous inhomogeneity apply:

⟨C:ε⟩_r = C_r : A_r ,      N_r = (C_r - C₀) : A_r .

true for every ellipsoid, cylinder and crack. false for internally heterogeneous inclusions such as LayeredSphere, whose average stress must be assembled layer by layer — no single C_r represents them, and feeding the phase property into the formulas above gives a wrong answer (it can even come out with the opposite sign).

Scheme kernels must branch on this trait rather than inlining the homogeneous formula.

is_homogeneous_inclusion(::LayeredSphere) -> false

A composite sphere has no single representative property: its average stress must be summed over layers. See stress_strain_loc.

is_homogeneous_inclusion(::LayeredSpheroid) -> false

A composite spheroid has no single representative conductivity: its average flux must be obtained from the transfer-matrix recurrence, not from (K₁ - K₀)·A. See flux_gradient_loc.

loc_and_stiffness(incl, P₁, P₀; kw...) -> (A, N)

Bundled evaluation of the dilute concentration tensor A (strain_strain_loc / gradient_gradient_loc) and the size-independent contribution tensor N (stiffness_contribution / conductivity_contribution) of the same inclusion in the same reference medium, sharing the single expensive Hill / recurrence solve.

Both quantities are needed per phase by Mori-Tanaka and by the self-consistent kernels; computing them separately evaluates hill_tensor — the dominant cost — twice with byte-identical arguments.

Internal seam — not exported.

loc_and_stress_average(incl, P₁, P₀; kw...) -> (A, B)

Bundled (strain_strain_loc, stress_strain_loc) — resp. (gradient_gradient_loc, flux_gradient_loc) — sharing one localization solve. Same bitwise contract as loc_and_stiffness.

Internal seam — not exported.

compliance_and_stiffness_contribution(incl, P₀; kw...) -> (H, N)

Bundled two-argument contribution pair of a flat inclusion — the density-amount counterpart of loc_and_stiffness, sharing whatever single expensive solve produces both (for a crack, one cod_tensor).

Returns (compliance_contribution, stiffness_contribution) for a 4th-order P₀ and (compliance_contribution, conductivity_contribution) for a 2nd-order one. Same bitwise contract as loc_and_stiffness: the generic fallback in src/contribution.jl is that pair, and any specialization is a pure performance concern.

Internal seam — not exported.

MeanFieldHom.CustomInclusions

The user-defined inclusion contract: everything needed to plug a morphology MeanFieldHom knows nothing about into every homogenization scheme, in elasticity as in transport.

Two exports:

• CustomInclusion — a concrete, callback-driven inclusion, the value-level analogue of subtyping AbstractCustomInclusion;
• check_inclusion_interface — a conformance checker that works on any AbstractInclusion, reporting which entry gate it satisfies and what is missing.

The contract itself (four levels, three entry gates, the amount × contribution seam) is specified in docs/src/developer/adding_inclusion.md; the user-facing tutorial is docs/src/manual/custom_inclusions.md.

Loaded after localization.jl and contribution.jl, because the fallbacks here invoke the generic methods defined there.

AbstractCustomInclusion{T} <: AbstractInclusion{T}

Supertype for user-supplied inclusion morphologies that do not fit any of the built-in families — non-ellipsoidal shapes, patterns whose response is obtained from an external solver (finite elements, series expansion, …), or hand-crafted approximate formulas. It is the MeanFieldHom counterpart of the user_inclusion extension point of the C++/Python echoes codebase.

Subtyping this abstract type is not mandatory — an inclusion may equally well subtype AbstractEllipsoidalInclusion, AbstractCrack or AbstractLayeredInclusion when it genuinely belongs to that family, and will then inherit that family's dispatch rules. What this branch buys is a neutral home for morphologies that belong to none of them: it is disjoint from the other three, so adding methods for it can never create a dispatch ambiguity, and hill_tensor accepts it directly (the ellipsoid-typed entry point does not).

See the developer page Adding a new inclusion for the full interface contract, and CustomInclusion for a ready-made concrete type driven by callbacks.

CustomInclusion{dim,T,B,NT} <: AbstractCustomInclusion{T}

Inclusion whose mechanical / transport response is supplied by user callbacks, so that an arbitrary morphology becomes a first-class citizen of every homogenize scheme without defining a new type.

Construction

CustomInclusion(; dim, semi_axes, basis, euler_angles, homogeneous,
                density_factor, <callbacks>...)
CustomInclusion(semi_axes; kw...)          # convenience

Callbacks — pick one entry gate

Every callback is a keyword argument named exactly after the generic it implements, and receives the tensor arguments only (the geometry is captured by the closure). All callbacks must accept trailing keyword arguments: the schemes forward method, tolerances, and — for density-based phases — K_interface / α_interface.

Gate A yields all eight localization and all four contribution tensors for free; gate B yields the derived localizations and the contributions. P₀ is a 4th-order stiffness in elasticity and a 2nd-order conductivity in transport — a single callback may serve both by dispatching on its argument.

Orientation averaging (IsoSymmetrize / TISymmetrize) is applied by the scheme after the callback returns, so the returned tensors must be expressed in the global frame.

Example — a sphere, through gate A

sphere = CustomInclusion((1.0, 1.0, 1.0);
    hill_tensor = (C₀; kw...) -> hill_tensor(Ellipsoid(1.0, 1.0, 1.0), C₀; kw...))

rve = RVE(:m, Dict(:C => iso_stiffness(10.0, 6.0)))
add_phase!(rve, :i, sphere, Dict(:C => iso_stiffness(1.0, 0.5)); fraction = 0.2)
homogenize(rve, MoriTanaka(), :C)

See also check_inclusion_interface and the developer page Adding a new inclusion for the full contract.

CustomShape

Neutral shape_trait tag carried by CustomInclusion. No kernel dispatches on it — a custom inclusion supplies its own response — but the interface requires some shape trait.

check_inclusion_interface(incl; physics = :elasticity, amount = :fraction,
                          verbose = true) -> Bool

Report which parts of the inclusion interface incl satisfies, and return true when it can be fed to the homogenization schemes for the requested physics (:elasticity or :conduction) and amount (:fraction or :density).

Checks, in order:

1. Level 0 — dimension, inclusion_basis, shape_trait, shape_tensor.
2. Level 1 — which entry gate is available: the Hill tensor, the localization tensor, or the contribution tensors.
3. Level 2 — for amount = :density, the three-argument delta_* seams.

Works on any AbstractInclusion, not only CustomInclusion — use it on your own type before wiring it into an RVE.

MeanFieldHom.FiniteElements

Inclusions whose response is obtained from a finite-element resolution of the Eshelby problem rather than from a closed form — the package's own demonstration that the CustomInclusions contract is enough to reach every scheme.

Both use the finite Eshelby cell with a first-order corrected boundary condition of Adessina, Barthélémy, Lavergne & Ben Fraj, Int. J. Eng. Sci. 119 (2017) 1-15: the infinite matrix is truncated to a ball of finite radius and the truncation bias is removed by adding the inclusion's own dipole far field to the imposed boundary displacement.

The types, the Fourier operators, the boundary data and the algebra of the corrected boundary condition all live here. What a package extension supplies is only the discretization — a mesh, scalar Lagrange spaces, an assembly and a quadrature — through the nine generics of FEBackend.

Two backends exist, and both serve both morphologies: FerriteBackend (import Ferrite, FerriteGmsh, Gmsh) and GridapBackend (import Gridap, GridapGmsh). An inclusion built without naming one takes AutoBackend and picks at its first solve; with neither loaded, that solve errors informatively.

See docs/src/manual/fe_inclusions.md and docs/src/applications/recycled_aggregate.md.

FECache()

Mutable side-store of a finite-element inclusion: the assembled discretization (built once, on first use) and the response tensors already computed, keyed on the reference medium.

assemblies counts the factorizations actually performed — used by the tests to prove the memoization works. The field is deliberately kept out of the struct's numeric fields so that it never interferes with ForwardDiff reconstruction of a geometry parameter.

fe_assembly_count(incl) -> Int

Number of finite-element assemblies (equivalently, factorizations) actually performed for incl so far. Every distinct reference medium costs one; a repeat costs none. Useful to check that the memoization of FECache is doing its job.

fe_reset!(incl) -> incl

Drop the cached discretization and every memoized response tensor.

FEBackend

Which finite-element library performs the solve. Concrete singletons: FerriteBackend, GridapBackend and the default AutoBackend.

The three are always defined, whether or not the corresponding package is loaded — naming a backend costs nothing, only solving with it requires the matching extension.

AutoBackend()

Default backend of a finite-element inclusion: pick whichever backend is loaded, at the first solve rather than at construction.

Deferring the choice is what lets an inclusion be built, printed, stored in an RVE and passed around in a session where no finite-element package has been imported; the informative error arrives only when a scheme actually asks for a localization tensor.

Priority is FerriteBackend then GridapBackend. Loading both is not an error — it just means the first is chosen. To pick deliberately, pass backend = GridapBackend() to the constructor.

FerriteBackend()

Solve with Ferrite.jl; needs import Ferrite, FerriteGmsh, Gmsh. The reference implementation, and the faster of the two to run.

GridapBackend()

Solve with Gridap.jl; needs import Gridap, GridapGmsh (GridapGmsh carries its own gmsh, so Gmsh.jl is not required on this path).

Gridap states the weak form directly — ∫( ε(v) ⊙ (σ∘ε(u)) )dΩ for the crack, ∫( Bᵐ(v)' * D * Bᵐ(u) * ρ )dΩ for the axisymmetric modes — which makes it the easier of the two to read and to modify.

_build_gmsh_axi_model(gmsh, a, a_core, d, R, h_in, h_out)

Populate the current gmsh session with the meridian half-plane of a sphere of radius a holding a core of radius a_core centred at z = d, itself embedded in a ball of matrix of radius R.

Element size is h_in on the inclusion and its core and h_out on the outer boundary, gmsh interpolating in between.

The gmsh module is passed in rather than imported: gmsh is a weak dependency, and the two finite-element backends reach it by different routes — Gmsh.gmsh for Ferrite, GridapGmsh.gmsh for Gridap, both over the same gmsh_jll. Taking it as an argument keeps this geometry here in src/, shared, instead of duplicating it in each extension.

The caller owns gmsh.initialize() / gmsh.finalize().

fe_axi_grid(backend, incl)

Backend-native mesh of the meridian half-plane of incl, carrying the cell sets "core", "shell", "matrix" and the boundary sets "outer", "axis" of _build_gmsh_axi_model. The first node coordinate is the cylindrical radius ρ, the second the axial coordinate z.

fe_axi_grid_counts(backend, grid) -> (; ncells, nnodes, ncells_by_set)

Cell and node counts, ncells_by_set being a Dict{String,Int} over the three regions. Diagnostics only.

fe_axi_region_volume(backend, grid, set) -> Float64

Volume of revolution 2π ∫_set ρ dρ dz of one region, on the geometric interpolation of the mesh. Diagnostics only — the driver measures its own volume with the mode's quadrature, and the two need not agree exactly.

fe_axi_mode(backend, grid, order, ncomp, axis_zeros) -> mode

Discretize one Fourier mode: ncomp scalar Lagrange fields of degree order on grid, with a quadrature exact to degree 2 * order + 1.

The dof numbering must span the whole space, with no Dirichlet elimination — the driver does the free/prescribed split itself, so that one factorization serves every right-hand side.

axis_zeros lists the component indices that the axis regularity of this mode forces to vanish on the set "axis".

fe_axi_dof_split(backend, mode) -> (ndofs, free, presc)

Total dof count and the two index vectors, in the numbering of fe_axi_mode. presc is the sorted union of the outer-boundary dofs and the axis-pinned dofs; free is its complement.

fe_axi_set_dirichlet!(backend, mode, u, f) -> u

Write the Dirichlet data of one right-hand side into u: u[d] = f(ρ_d, z_d)[k] for every dof d of the outer boundary, k being its component, then u[d] = 0 for every axis-pinned dof.

Two things make this the delicate function of the contract.

The order matters. The poles (0, ±R) belong to both the "outer" and the "axis" sets. The axis must win, so the zeroing comes second.

It must not touch the matrix. The driver assembles and factorizes once, then calls this once per load case. Re-assembling here would multiply the cost of a solve by the number of loads.

fe_axi_stiffness(backend, mode, Dmap, Bop) -> AbstractMatrix

Stiffness of one mode over the whole dof numbering,

K = Σ_regions ∫_region Bop(v)' * D_region * Bop(u) * ρ dρ dz .

Dmap is a Vector{Pair{String,Matrix{Float64}}} — a vector, not a Dict, so that the assembly order is reproducible — mapping a cell-set name to that region's material matrix in the cylindrical (ρ, θ, z) basis.

Bop(N, dNρ, dNz, ρ) -> Matrix{Float64} of size nrow × ncomp is the generalized-strain operator of one scalar shape function, already closed over the Fourier mode. It is R-linear in (N, dNρ, dNz), so a backend that manipulates whole trial functions rather than shape functions may apply it to (u_c, ∂ρu_c, ∂zu_c) directly instead of building an element B matrix.

The ρ in the measure is the single factor that turns a plane problem into a solid of revolution.

fe_axi_average(backend, mode, Dmap, u, Bop, proj, sets) -> (prim, dual, V)

Volume averages over ∪ sets of the generalized strain and of the associated generalized stress, both projected by proj onto the Kelvin basis of the mode, plus the volume of revolution V = 2π ∫ ρ dρ dz of that union:

prim = ∫ proj(B u) ρ / ∫ ρ ,      dual = ∫ proj(D · B u) ρ / ∫ ρ .

The azimuthal integration has already been performed analytically inside proj; what remains is the meridian quadrature.

_resolve_backend(b) -> FEBackend

Turn AutoBackend into a concrete backend, and check that a concrete one is actually available. Called once per inclusion, at the first solve; the result is pinned in the cache, so a single inclusion never mixes two backends' grids.

_build_gmsh_crack_model(gmsh, a, b, R, htipdiv)

Populate the current gmsh session with the crack-in-a-ball model: a ball of radius R centred on an elliptical crack of semi-axes (a, b) lying in the z = 0 plane, refined to min(a,b)/htipdiv in a torus hugging the crack front and coarsening to R/3 at the outer boundary.

The gmsh module is passed in rather than imported, for the reason given in _build_gmsh_axi_model: it is a weak dependency reached by a different route from each backend. The caller owns gmsh.initialize() / gmsh.finalize().

_weld_msh_crack_front(path, a, b; ell_tol, z_tol) -> Int

Merge the duplicated nodes of the crack front in a written MSH 4.1 file, in place, and return how many pairs were welded.

The Crack plugin duplicates the nodes of the front along with those of the lips, in spite of OpenBoundaryPhysicalGroup (still true in gmsh 4.15). Left alone the crack is effectively half an element longer than asked for, and the opening comes out 10-20 % too large. The lips must stay split — that discontinuity is the crack — so a blanket removeDuplicateNodes is not an option: only nodes on the ellipse (x/a)² + (y/b)² = 1, z = 0 are merged.

Working on the file rather than on the live gmsh model is what makes this shared: a node merge is a renumbering of the element connectivity, and every backend reads the same file. Nodes left unreferenced are harmless — both mesh readers drop them.

fe_crack_grid(backend, crack)

Backend-native mesh of the ball holding the crack, with the boundary sets "outer" (the sphere) and "crack" (both lips, whose nodes the gmsh Crack plugin has duplicated). The crack front must already be welded — see _weld_msh_crack_front.

fe_crack_counts(backend, grid) -> (; ncells, nnodes, nfacets_up, nfacets_dn,
                                     area_up, area_dn)

Mesh diagnostics: cell and node counts, and the facet count and area of each lip. The two areas must both equal πab — that is what says the plugin split the surface cleanly and the front weld did not glue the lips back together.

fe_crack_space(backend, grid, order) -> space

A vector-valued Lagrange space of degree order over the whole mesh, with a quadrature exact to degree 2 * order and no Dirichlet elimination — the driver splits the dofs itself so that one factorization serves all six right-hand sides.

The lips need no special treatment: they are traction-free naturally, because their nodes are distinct. There is no interface term, no multiplier and no contact condition anywhere in this problem.

fe_crack_dof_split(backend, space) -> (ndofs, free, presc)

Total dof count and the two index vectors, presc being the dofs of the outer sphere.

fe_crack_set_dirichlet!(backend, space, u, f) -> u

Write u[d] = f(x_d)[k] for every dof d of the outer sphere, x_d its node and k its component. Called once per right-hand side; must not touch the matrix.

fe_crack_stiffness(backend, space, C) -> AbstractMatrix

Stiffness of linear elasticity, ∫ ε(v) : ℂ : ε(u) dΩ, over the whole dof numbering.

C is a Tensors.SymmetricTensor{4,3} and is isotropic: the corrected boundary condition uses the closed-form Kelvin dipole field, so the driver refuses anything else long before reaching here. A backend may therefore work from (λ, μ) instead of from the full tensor.

fe_crack_mean_jump(backend, space, u, S_f, b) -> Vector{3}

⟨[[u]]⟩ / b, the opening averaged over the crack surface and normalized by the semi-minor axis — the convention of cod_tensor.

Measured as a surface integral of the trace of u on each lip, with no assumption on the opening profile. The two lips are told apart by the sign of n ⋅ e₃, n being the outward normal of the adjacent element: the lip whose element sits above the crack carries n = -e₃ and contributes +u.

FEEllipticCrack(a, b; euler_angles = (), radius_ratio = 5.0,
                htipdiv = 12.0, order = 2)

Flat elliptical crack whose crack-opening-displacement tensor is computed by finite elements instead of the closed form of EllipticCrack.

It subtypes AbstractCrack and declares the standard shape_trait, so implementing cod_tensor is all it takes: ℍ, ℕ, 𝐑, 𝐍K, the bundled pair and the four `delta*with the Budiansky4π/3prefactor are inherited. It is a drop-in replacement forEllipticCrack` in every scheme — the point of the exercise being that the same machinery accepts a morphology for which no closed form exists.

Requires Ferrite, FerriteGmsh and Gmsh to be loaded.

Scope

Isotropic matrix only. An anisotropic reference medium needs the anisotropic Green-function gradient (Pan-Chou, or the Barnett-Willis line integral), which is not implemented — see docs/src/developer/roadmap.md. Isotropy is tested on the tensor's content, not its TensND type.

That restriction bites in the iterative schemes: SelfConsistent and DifferentialScheme re-evaluate the crack in the current estimate, which for a family of parallel cracks is transversely isotropic. Add symmetrize = IsoSymmetrize() to the phase and the scheme hands the kernel an isotropic reference at every iteration — which also collapses the whole orientation family onto one cached solve.

ForwardDiff cannot be propagated through the solve (the sparse factorization is Float64-only), so use finite differences for sensitivities.

Example

using MeanFieldHom, Ferrite, FerriteGmsh, Gmsh

crack = FEEllipticCrack(1.0, 0.25; htipdiv = 12.0)
C₀ = iso_stiffness(0.8333, 0.3846)          # E = 1, ν = 0.3

B_fe = cod_tensor(crack, C₀)                 # finite elements
B_an = cod_tensor(EllipticCrack(1.0, 0.25), C₀)   # closed form

rve = RVE(:M)
add_matrix!(rve, Ellipsoid(1.0), Dict(:C => C₀))
add_phase!(rve, :cracks, crack, Dict(:C => C₀); density = 0.05,
           symmetrize = IsoSymmetrize())
homogenize(rve, MoriTanaka(), :C)

See also FEMeshOptions, FECache.

FEMeshOptions(; radius_ratio = 5.0, htipdiv = 12.0, order = 2)

Discretization settings of a finite-element inclusion.

order = 2 is strongly recommended: linear tetrahedra badly under-resolve the square-root field at the crack front.

fe_mesh_report(crack) -> NamedTuple

Mesh diagnostics: cell, node and dof counts, the two lip facet counts and their measured areas against the exact πab. Builds the discretization if it does not exist yet, and caches it.

Both lip areas equalling πab is what says the Crack plugin split the surface cleanly and the front weld did not glue the lips back together.

fe_cod_breakdown(crack, C₀) -> (; B_s, B_u, B_inf, B_s_glob, B_inf_glob)

Diagnostic view of the corrected solve: the COD tensor of the finite cell B_s, the response B_u to the crack's own dipole far field, and the infinite-medium result B_inf = (1 - B_u)⁻¹ B_s, all in the crack's local frame, plus B_s and B_inf rotated back to the global frame.

norm(B_u) measures how much work the boundary correction is doing; it should fall like (a/R)³, and B_inf — unlike B_s — should be insensitive to radius_ratio. That contrast is the practical proof that the correction is wired correctly.

Bypasses the cache.

FEExcenteredSphere(a, (P_core, P_shell); core_fraction, eccentricity = 0.0,
                   euler_angles = (), radius_ratio = 4.0, nradial = 24,
                   coarsening = 6.0, order = 2)

Spherical inclusion of radius a containing a spherical core placed off the centre, resolved by axisymmetric Fourier finite elements.

Geometry

α = 0 is the concentric two-layer sphere, for which LayeredSphere gives the exact Hervé-Zaoui answer — the reference this type is validated against. α → 1 brings the core tangent to the outer surface. The symmetry axis is axis; the response is transversely isotropic about it (and isotropic at α = 0).

What it provides

The inclusion is heterogeneous, so it enters through gate B of the inclusion contract with both localization tensors — the strain-side A_εε and the stress-side A_σε (resp. A_∇∇ and A_q∇ in transport). All schemes that consume those two follow: Dilute, MoriTanaka, Maxwell, PonteCastanedaWillis, SelfConsistent, DifferentialScheme.

Voigt and Reuss work too, which is not automatic for a heterogeneous inclusion: a bound averages the constituent properties over the inclusion and therefore needs its internal volume fractions, which the RVE does not carry. Here the geometry fixes them exactly — w for the core, 1 - w for the shell, whatever the eccentricity — so the type implements Schemes._layer_voigt / Schemes._layer_reuss and the bounds follow.

Properties

Like LayeredSphere, the constituent properties live in the geometry object, core first then shell, and the Dict handed to add_phase! is a placeholder that the kernel ignores. Build one object per physics: a pair of Tens{4,3} for elasticity, a pair of Tens{2,3} for conduction.

Requires a finite-element backend to be loaded — Ferrite, FerriteGmsh and Gmsh, or Gridap and GridapGmsh; pass backend = GridapBackend() to pick the second when both are available, and see FEBackend. The reference medium must be isotropic (the corrected boundary condition uses the closed-form Kelvin dipole field); the constituents may be isotropic or transversely isotropic about the symmetry axis.

Example

using MeanFieldHom
import Ferrite, FerriteGmsh, Gmsh        # or: import Gridap, GridapGmsh

C_core, C_shell = iso_stiffness(20.0, 12.0), iso_stiffness(6.0, 4.0)
C₀ = iso_stiffness(10.0, 6.0)
incl = FEExcenteredSphere(1.0, (C_core, C_shell);
                          core_fraction = 0.5, eccentricity = 0.4)

A, B = fe_axi_localization(incl, C₀)          # both tensors, one solve

rve = RVE(:m, Dict(:C => C₀))
add_phase!(rve, :rca, incl, Dict(:C => C₀); fraction = 0.3)
homogenize(rve, MoriTanaka(), :C)

See also FEAxiMeshOptions, fe_axi_breakdown, fe_axi_mesh_report.

FEAxiMeshOptions(; radius_ratio = 4.0, nradial = 24, order = 2,
                 coarsening = 6.0)

Discretization settings of an axisymmetric finite-element inclusion.

The mesh is two-dimensional (the meridian half-plane), so refining is cheap: nradial = 40 on a triangle mesh still solves in a fraction of a second.

fe_axi_localization(incl, P₀) -> (A, B)

The pair of localization tensors of an axisymmetric finite-element inclusion, computed in one solve — the strain-side A_εε and stress-side A_σε in elasticity, A_∇∇ and A_q∇ in transport.

Calling the two generics separately returns exactly the same tensors at no extra cost: they share the memoized solve.

fe_axi_breakdown(incl, P₀) -> NamedTuple

Diagnostic view of the corrected axisymmetric solve. Besides the corrected localization tensors A, B it returns their uncorrected counterparts — those of the plain truncated cell, u|∂Ω = E·x — and the per-mode blocks (A_E, B_E, A_p, B_p, A, B, X), plus the measured inclusion volume.

A_uncorrected drifts with radius_ratio like (a/R)³ while A does not: that contrast is the practical proof that the correction is wired correctly. Bypasses the cache.

fe_axi_mesh_report(incl) -> NamedTuple

Mesh diagnostics: cell and node counts per region, and the measured volumes of the core, the shell and the whole cell against their exact values. Builds the grid if it does not exist yet, and caches it.

The volumes are measured on the geometric interpolation with a quadrature of its own, so they do not match the solver's volume to the last digit; both converge to the same limit.

core_radius(incl) -> Real

Radius of the off-centre core, a·w^(1/3).

core_offset(incl) -> Real

Distance from the inclusion centre to the core centre, α·(a − a_core).

tensor_order(incl) -> Int

4 for an elasticity object (constituents are stiffness tensors), 2 for a transport one (conductivity tensors). Fixed at construction by the pair of properties handed in.

ExcenteredSphereShape

shape_trait of FEExcenteredSphere. No kernel dispatches on it — the inclusion supplies its own localization tensors.

green_gradient_iso(C₀::TensISO{4,3}, x) -> SArray{Tuple{3,3,3}}

Gradient $\partial G_{ij}/\partial x_k$ of the Kelvin Green function of an isotropic elastic matrix C₀, evaluated at x ≠ 0.

With $r = \|x\|$, $\hat n = x/r$ and $A = 1/\bigl(16\pi\mu(1-\nu)\bigr)$ the Kelvin solution reads

\[G_{ij}(x) = \frac{A}{r}\left[(3-4\nu)\,\delta_{ij} + n_i n_j\right],\]

so that

\[\frac{\partial G_{ij}}{\partial x_k} = \frac{A}{r^{2}}\left[-(3-4\nu)\,\delta_{ij} n_k + \delta_{ik} n_j + \delta_{jk} n_i - 3\,n_i n_j n_k\right].\]

Returned as a static 3×3×3 array indexed [i, j, k].

Throws a DomainError at the origin. Type-generic (Float64, ForwardDiff.Dual, symbolic scalars).

See also dipole_displacement_iso.

dipole_displacement_iso(C₀::TensISO{4,3}, x, Π) -> SVector{3}

Displacement field at x generated in an infinite isotropic matrix C₀ by a point polarization (force dipole) Π, i.e. the contraction

\[u_i(x) = \frac{\partial G_{ij}}{\partial x_k}(x)\; \Pi_{jk}.\]

Π has the dimension of a stress times a volume: for an inclusion of volume $V_{\mathcal I}$ carrying a uniform polarization $\pi$ (i.e. $\sigma = \mathbb C_0 : \varepsilon + \pi$ inside it), $\Pi = V_{\mathcal I}\,\pi$.

For a symmetric Π the closed form collapses to

\[u(x) = \frac{1}{16\pi\mu(1-\nu)r^{2}} \Bigl[-2(1-2\nu)\,\Pi\!\cdot\!\hat n + \mathrm{tr}(\Pi)\,\hat n - 3(\hat n\!\cdot\!\Pi\!\cdot\!\hat n)\,\hat n\Bigr],\]

which is the form evaluated here — it costs one matrix-vector product instead of building the full 3×3×3 gradient, and is the expression used for the corrected boundary condition of a finite Eshelby cell.

Π may be given as a 3×3 matrix or as a TensND 2nd-order tensor; in the latter case it is read in the canonical frame.

See also green_gradient_iso.

MeanFieldHom.NeuralInclusions

Inclusions whose response is produced by a trained neural network instead of a closed form or a finite-element solve — a fourth route into the CustomInclusions contract, alongside the analytic families, the layered patterns and FiniteElements.

Why bother, when the analytic Hill tensor is exact

For the ellipsoid the surrogate is not faster than the closed form, and it is less accurate. What it buys is two things the expensive routes cannot give:

• differentiability. A surrogate is a smooth function of its inputs, so derivative(rve, scheme, geometry(:phase, :field)) reaches a morphology parameter. The finite-element inclusions refuse that request outright (FiniteElements.jl): their solve runs in Float64 and memoizes on the reference medium, so the derivative would come out as a silent zero.
• cost, once the teacher is expensive. One FEExcenteredSphere evaluation is three assemblies and eight solves, and an iterative scheme changes the reference medium at every iteration, defeating the cache. A surrogate trained on that solve answers in microseconds.

The ellipsoid is therefore the validation case, not the application: it is the one morphology where the labels are exact, so every part of the pipeline can be checked against a closed form before being pointed at something unknown.

The pieces

• mlp.jl — a dependency-free, type-generic multilayer perceptron. Evaluation needs LinearAlgebra and nothing else, which is what lets a committed model be loaded anywhere and be traversed by ForwardDiff.
• specs.jl — the physics: which tensor class is predicted, and how the exact invariances (symmetry class, major symmetry, homogeneity in the reference moduli, frame indifference) are enforced rather than fitted.
• surrogate.jl — NeuralSurrogate: weights, standardization, output specification, validity box and Provenance.
• dataset.jl — Halton sampling of a SampleBox and labelling by a teacher, the one seam that changes between morphologies.
• training.jl — TrainingOptions and the fallback train_surrogate; the optimizer lives in MeanFieldHomLuxExt.
• io.jl — save_surrogate / load_surrogate, JSON.

Fitting needs Lux, Zygote and Optimisers (weak dependencies, as for the finite-element backends); evaluating needs none of them.

See docs/src/manual/neural_inclusions.md and scripts/84_neural_inclusion_ellipsoid.jl.

NeuralHillInclusion(semi_axes; elastic = nothing, transport = nothing,
                    basis = nothing, euler_angles = (), guard = :warn)

Inclusion whose Hill tensor is produced by a trained NeuralSurrogate — entry gate A of the inclusion contract, so all eight localization tensors, all four contribution tensors and every homogenization scheme follow.

At least one surrogate is required, and each must match its tensor order. A surrogate is checked against the geometry at construction: a HillTI model needs a spheroid, a HillISO one a sphere, so a mismatch is a constructor error rather than a wrongly shaped tensor much later.

Example

s = load_surrogate(NeuralInclusions.model_path("spheroid_hill_iso_elastic"))
incl = NeuralHillInclusion((1.0, 1.0, 0.2); elastic = s)

rve = RVE(:m, Dict(:C => iso_stiffness(20.0, 12.0)))
add_phase!(rve, :i, incl, Dict(:C => iso_stiffness(60.0, 30.0)); fraction = 0.2)
homogenize(rve, MoriTanaka(), :C)

See also NeuralLocalizationInclusion, and docs/src/manual/neural_inclusions.md.

NeuralLocalizationInclusion(semi_axes; strain = nothing, stress = nothing,
                            gradient = nothing, flux = nothing,
                            fractions = nothing, properties = nothing,
                            basis = nothing, euler_angles = (), guard = :warn)

Inclusion whose localization tensors are produced by trained surrogates — entry gate B, the only way in for an internally heterogeneous morphology, which has no Hill tensor at all. This is the shape a surrogate trained on fe_axi_localization takes.

Because is_homogeneous_inclusion is false, gate B costs two tensors per physics: the strain side and the stress side, since 𝔸_σε = ℂ₁:𝔸_εε presupposes a single uniform ℂ₁ that a heterogeneous inclusion does not have. Given both, the generic contributions switch to the exact ℕ = 𝔸_σε − ℂ₀:𝔸_εε and ℍ = (𝔸_εε − 𝕊₀:𝔸_σε):𝕊₀, so gate B is a complete entry point.

A gate-B surrogate's features are the morphology parameters and the contrast ratios — never absolute moduli. The reason is that the constituents live inside the inclusion, so scaling ℂ₀ alone changes the contrast: the exact invariance is under a simultaneous scaling of the reference medium and of every constituent, which leaves 𝔸_εε unchanged and multiplies 𝔸_σε by the factor. The ℙ(λℂ₀) = ℙ(ℂ₀)/λ homogeneity that gate A exploits does not transfer.

Supplying only one tensor of a pair is refused at construction: the omission is silent otherwise — Dilute and MoriTanaka stay right while SelfConsistent and AsymmetricSelfConsistent, the schemes that consume the stress average, drift by several percent.

NeuralShape

Neutral shape_trait of the neural inclusions. No kernel dispatches on it — a surrogate supplies its own response — but the contract requires some shape trait.

NeuralSurrogate

A trained network plus its physical contract: which features it consumes, how both ends are standardized, which tensor it produces, and over what box it is entitled to be believed.

Construction

Built by train_surrogate! or read back by load_surrogate. The direct constructor is keyword-only and validates every length against the network's own input and output widths, because a surrogate whose feature list and input layer disagree is a silent mis-evaluation rather than an error.

Evaluation

(s::NeuralSurrogate)(x_raw, P₀, frame; guard = :warn)

x_raw is the unstandardized feature vector — the inclusion builds it, since only the inclusion knows its own geometry (see _raw_features). frame is the symmetry axis or material frame passed straight to build.

Generic in the element type of x_raw: a ForwardDiff.Dual feature yields a Dual tensor, which is what makes a morphology sensitivity possible.

See also Provenance, validate_surrogate, save_surrogate.

Provenance

Where a surrogate came from and how well it does, recorded at training time and serialized with the weights.

max_block_error is the number a test tolerance should be derived from — never a hard-coded literal, so that a retraining cannot silently loosen a threshold. The per-component vectors are diagnostics: a structurally vanishing component (𝕍ᴬ has no ℓ₃) has no meaningful relative error of its own, which is exactly why the headline number is measured against the block.

worst_error(p::Provenance) -> Float64

Worst error recorded on the held-out set, relative to the tensor's own magnitude (max_block_error), or Inf when the surrogate has never been validated — so that a tolerance derived from it fails loudly rather than passing by accident.

check_domain(s, x_raw, guard) -> Nothing

Verify that every feature lies inside the box the surrogate was trained on.

guard is :error (refuse), :warn (warn once per call site, the default) or :none (trust the caller). Extrapolating a network is not a graceful degradation — the error can be arbitrary and there is no diagnostic in the result — so the default is deliberately noisy.

The bounds are inclusive up to a relative slack of 1e-9 of the box width, so that reconstructing a bound by a slightly different arithmetic path — log(0.1) against a stored log(0.10) — does not trip the guard on the last bit.

predict_components(s, x_raw; guard = :none) -> Vector

The surrogate's raw prediction — the dimensionless or shape components, before any contraction with the reference moduli. Used by validate_surrogate, which compares against labels living in exactly this space, and by training diagnostics.

AbstractHillClass

Which structured TensND type a surrogate predicts, and therefore how many components it emits. A class fixes the tensor order too: order 4 is elasticity, order 2 is transport.

Isotropic 4th-order tensor α𝕁 + β𝕂 — the Hill tensor of a sphere.

Transversely isotropic, major-symmetric 4th-order tensor: the five Walpole components (ℓ₁, ℓ₂, ℓ₃, ℓ₅, ℓ₆) of TensTI{4,·,5} — the Hill tensor of a spheroid in an isotropic matrix.

Orthotropic 4th-order tensor: the nine components (C₁₁, C₂₂, C₃₃, C₁₂, C₁₃, C₂₃, C₄₄, C₅₅, C₆₆) of TensOrtho — the Hill tensor of a triaxial ellipsoid in an isotropic matrix.

Isotropic 2nd-order tensor — the transport Hill tensor of a sphere.

Transversely isotropic 2nd-order tensor (a, b) of TensTI{2,·,2} — the transport Hill tensor of a spheroid.

AbstractOutputSpec

How the network's raw output vector becomes a Hill tensor. Two implementations, differing only in how much of the material dependence is exact rather than learned: DimensionlessHill and AffineHill.

DimensionlessHill(class)

The network predicts the ncomponents(class) components of the dimensionless Hill tensor scale · ℙ. Exact in the scale of the reference moduli and in the symmetry class; ν₀ is an input feature and its dependence is learned.

The general-purpose choice: it needs nothing of the reference medium beyond a scale, so the same shape of surrogate transfers to an anisotropic matrix, and to a localization pair (gate B) where no affine structure exists.

AffineHill(class)

The network predicts the components of the shape-only tensors 𝕌ᴬ and 𝕍ᴬ — nterms(class) · ncomponents(class) numbers — and the decoder contracts them with material_coeffs. The whole material dependence is then exact: ν₀ is not an input, and the surrogate is a function of the shape alone.

Only available for an isotropic reference medium, which is where the affine structure comes from.

ncomponents(class) -> Int

Number of independent components the class carries — the width of one term of the network output.

tensor_order(class) -> Int

4 for elasticity, 2 for transport. The localization and contribution generics are declared per order, so this is what decides which physics a surrogate serves.

tensor_order(t::AbstractTens) -> Int

Order of a TensND tensor, read off its type — the same notion as for a class, so the two can be compared directly.

nterms(spec) -> Int

How many shape tensors the network predicts per component: 1 for DimensionlessHill, and for AffineHill the length of material_coeffs — 2 in elasticity, 1 in transport.

noutputs(spec) -> Int

Width of the network's output layer.

needs_nu(spec) -> Bool

Whether the feature set has to carry ν₀. False for AffineHill, whose material dependence is exact, and for any transport surrogate.

build(class, c, frame) -> AbstractTens

Assemble the structured tensor of class from the component vector c. frame is the symmetry axis (an NTuple{3}) for the TI classes, the material frame (a TensND basis) for HillOrtho, and ignored for the isotropic ones.

components(class, P, frame; atol = 1.0e-8) -> NTuple

The independent components of P, in the order build expects — the exact inverse of build.

Goes through TensND.proj_tens rather than reading get_data off the tensor, for two reasons. First, robustness: a teacher is free to return whichever concrete type it likes, and the conduction kernels do exactly that — the analytic transport Hill tensor of a spheroid comes back as a generic Tens{2,3}, not a TensTI{2,·,2}, so reading its data would yield nine numbers where two were wanted. Second, and more valuable: proj_tens also returns the relative residual of the projection, which for a tensor genuinely in the class is zero to round-off. Checking it against atol turns a wrong axis, a wrong frame or a wrong class — all of which would otherwise train happily on corrupted labels — into a loud error at dataset-generation time.

decode(spec, z, P₀, frame) -> AbstractTens

Turn the network's untransformed output z into the Hill tensor, in the frame given.

z has length noutputs; for AffineHill it is read as a ncomponents × nterms column-major block, one column per shape tensor.

material_coeffs(class, P₀) -> Tuple

Coefficients of the exact affine decomposition of the Hill tensor on the shape-only tensors, for an isotropic reference medium:

• order 4: (d, 1/μ₀) with d = 1/(λ₀+2μ₀) − 1/μ₀, so that ℙ = d·𝕌ᴬ + (1/μ₀)·𝕍ᴬ;
• order 2: (1/k₀,), so that ℙ_K = 𝕍ᴬ/k₀.

Used by AffineHill. The number of entries is the number of terms the network has to predict per component.

dimensionless_scale(class, P₀) -> Number

The modulus by which ℙ is multiplied to make it dimensionless — 2μ₀ in elasticity, k₀ in transport. Used by DimensionlessHill: because ℙ is homogeneous of degree −1 in the reference moduli, scale · ℙ depends on the shape and on ν₀ alone, and on nothing at all in transport.

hill_class(name::Symbol) -> AbstractHillClass

Class from its serialized name: :iso, :ti, :ortho, :iso2, :ti2.

output_spec(name::Symbol, class::Symbol) -> AbstractOutputSpec

Rebuild an output specification from its serialized names.

apply_transform(kind, z)

Map a physical component to the space the network is fitted in (:identity or :log).

invert_transform(kind, y)

Inverse of apply_transform — network space back to a physical component.

_feature(::Val{name}, incl, P₀) -> Number

One raw (unstandardized) feature of incl under reference medium P₀.

Three conventions here are load-bearing, and each of them is a bug if broken.

The logarithm is not cosmetic: an aspect ratio's interesting range spans decades, and ω and 1/ω are the same amount of anisotropy, which only the logarithm makes symmetric.

:log_aspect is measured on the distinct axis, not on a fixed slot. Semi-axes are stored sorted descending (as Ellipsoid does), so a prolate spheroid is (ω, 1, 1) and an oblate one (1, 1, ω): log(a₃/a₁) would be negative for both and conflate the two families. Distinct-over-equal is instead a bijection onto the whole real line, with the sphere at the origin.

The triaxial pair is (a₂/a₁, a₃/a₂), not (a₂/a₁, a₃/a₁). With the sorted convention a₁ ≥ a₂ ≥ a₃, the admissible set of the first pair is exactly the box log_r2 ≤ 0, log_r32 ≤ 0, whereas the second pair is confined to a triangular wedge that no SampleBox can express.

raw_features(incl, s::NeuralSurrogate, P₀) -> Vector

The feature vector s expects, read off incl and P₀. Type-generic: a ForwardDiff.Dual semi-axis yields a Dual feature vector, and the whole chain downstream follows.

_class_frame(class, geom) -> frame

The frame argument build and components need: the symmetry axis for a TI class, the material frame for an orthotropic one, nothing for an isotropic one.

For the TI classes the column carrying the symmetry axis is derived from the semi-axes — column 1 for a prolate spheroid, column 3 for an oblate one, which is what the analytic kernels of Elasticity._hill_3d_iso use.

_canonical_axes(axes, basis) -> (axes, basis)

Semi-axes sorted descending with the basis columns permuted to match — the very convention Ellipsoid applies, through the very same helper.

Sharing it is not tidiness, it is correctness: a surrogate is trained on the components the analytic teacher returns, and those are expressed in the sorted frame. An inclusion that stored (1, 1, 3) where the teacher saw (3, 1, 1) would feed the network a reciprocal aspect ratio and read the answer in the wrong frame.

A 2-tuple or any non-3 dimension is returned untouched: the shipped surrogates are three-dimensional, and there is no sorted convention to honour elsewhere.

_spheroid_axis_index(a) -> Int

Which of the three semi-axes is the distinct one — the symmetry axis of a spheroid. Refuses a shape that is not a spheroid, because the alternative is to silently pick an axis and rotate the tensor by 90°.

SampleBox(names, lo, hi; scale = :linear)

Axis-aligned sampling box over the raw feature space, one entry per feature.

scale is per feature: :linear samples uniformly between lo and hi, :log samples uniformly in the logarithm — the right choice for an aspect ratio, whose interesting behaviour is spread over decades rather than over an interval.

The box travels into the trained surrogate as its check_domain limits, which is why it is a first-class object and not just a pair of loops in a script.

Example

# a spheroid aspect ratio over two decades, plus the matrix Poisson ratio
SampleBox([:log_aspect, :nu0], [log(1/20), 0.0], [log(20), 0.49])

Note that lo/hi are given in feature units: :log_aspect is already a logarithm, so scale = :linear on it samples the exponent uniformly, which is what a log sweep of the aspect ratio means.

Dataset

Raw features and untransformed targets, one column per sample.

Z lives in the space the AbstractOutputSpec predicts — dimensionless components for DimensionlessHill, shape-tensor components for AffineHill — not in the space of the tensor itself. Training and validate_surrogate both work there, so the reference moduli never enter the loss.

generate_dataset(geometry, response, spec, box, n; nvalidation = 0)
    -> (train::Dataset, validation::Dataset)

Sample box, label every point, and return the training and held-out sets.

The teacher is split in two callbacks, and the split is what makes the frame convention safe:

• geometry(x_shape) -> geom builds the morphology from the shape features (the box's features minus :nu0). geom needs only to answer semi_axes and inclusion_basis — an Ellipsoid does, and so does any user type.
• response(geom, P₀) -> AbstractTens evaluates the tensor to be learned. For the pilot this is hill_tensor; for a heterogeneous morphology it is one of the localization tensors, and for an expensive one it is a solve.

Because the frame in which the components are read is then obtained from geom through the very same _class_frame the inclusion uses at evaluation time, a mismatch between the two is impossible by construction — rather than being a silent 90° rotation discovered much later.

geometry(x) = Ellipsoid(1.0, 1.0, exp(x[1]))
response(g, C₀) = hill_tensor(g, C₀)
train, val = generate_dataset(geometry, response, spec, box, 4000; nvalidation = 1000)

What the labels are depends on spec:

• DimensionlessHill — one response call per sample, at the reference medium built from the sampled ν₀ (elasticity) or at unit conductivity (transport); the target is scale · ℙ.
• AffineHill — two response calls per sample at two different Poisson ratios, and the shape tensors 𝕌ᴬ, 𝕍ᴬ are recovered by solving the exact 2×2 affine system componentwise. ν₀ must not be in the box: the whole point is that it is not a degree of freedom.

sample_box(box, n; offset = 0) -> Matrix{Float64}

length(box) × n matrix of raw feature vectors, from the Halton sequence.

offset skips that many points, which is how a held-out set is drawn: pass offset = n_train and the two sets are disjoint yet drawn from the same well-distributed sequence.

grid_box(box, npts) -> Matrix{Float64}

Full tensor grid with npts points per feature — length(box) × npts^d columns. Used for error maps and for hitting the corners of the box, which a low-discrepancy sequence does not.

halton(i, base) -> Float64

Radical-inverse of i in base — the i-th point of the one-dimensional van der Corput sequence, in (0, 1).

feature_index(box, name) -> Int

Position of a named feature in the box, with a readable error when absent.

fit_scaling(train::Dataset; log_threshold = 30.0)
    -> (; x_shift, x_scale, y_kind, y_shift, y_scale)

Standardization of both ends, and the per-component output transform, read off the training set only — the held-out set must not inform them.

A component is fitted on a :log scale when it is strictly positive throughout and its dynamic range max/min exceeds log_threshold. That is the case of the oblate Walpole components as the aspect ratio goes to zero, which span decades and would otherwise monopolize a mean-squared loss.

validate_surrogate(s, data::Dataset)
    -> (; max_rel_error, rms_rel_error, max_block_error, worst)

Error of s over data, in the space the network predicts, at two granularities.

max_block_error — the headline number. For each sample, the ∞-norm of the prediction error over all components, relative to the ∞-norm of the target:

max_j ‖ẑⱼ − zⱼ‖_∞ / ‖zⱼ‖_∞ .

This is the quantity a test tolerance should be derived from, because it is the one that propagates: the decoded tensor is a linear combination of the whole component vector, so what matters is the error relative to the tensor's own magnitude — not to that of its smallest entry. worst is an alias for it.

max_rel_error / rms_rel_error — the per-component diagnostic. Component i on sample j is scored as

|ẑᵢⱼ − zᵢⱼ| / max(|zᵢⱼ|, floorᵢ) .

The floor is the component's own RMS over the set, so a component that merely passes through zero — the Walpole ℓ₃ of a near-spherical inclusion — does not report an unbounded error. When a component is identically zero over the whole set the floor falls back to the global RMS instead, which turns a 0/0 into the honest statement "this component is that fraction of the tensor's magnitude". A structurally-vanishing component is not a rarity: 𝕍ᴬ has no ℓ₃ at all, since the analytic kernel gives p₃ = d·u₃ with no 1/μ₀ term.

report_surrogate([io], s, data::Dataset; labels = nothing)

Print the per-component validation table of s over data: the relative error in the space the network predicts, component by component, plus the worst entry.

Dependency-free, so it works wherever a surrogate can be evaluated. labels overrides the component names, which are otherwise 1:n.

component_labels(class) -> Vector{Symbol}

Human-readable names of a class's independent components, for report_surrogate.

component_labels(spec::AbstractOutputSpec) -> Vector{Symbol}

Names of the network's outputs. For AffineHill the class labels are repeated once per shape tensor, suffixed 𝕌 and 𝕍.

TrainingOptions(; kw...)

Everything the optimizer needs, with defaults sized for the ellipsoid pilot (a few thousand samples, a few thousand parameters, seconds of wall time).

train_surrogate(spec, box, train, validation; options = TrainingOptions(),
                teacher_name = "", notes = "", history = nothing)
    -> NeuralSurrogate

Fit a surrogate of output specification spec over the sampling box, on the datasets produced by generate_dataset, and return it with its Provenance filled in from the held-out set.

Requires the training extension. Run

import Lux, Optimisers, Zygote

before calling. Without them this fallback method raises: evaluation of an already-trained surrogate needs none of the three, so they are weak dependencies rather than dependencies (scripts/nn/ carries an environment that has them).

Pass a Vector as history to have the learning curve recorded into it, one (; epoch, train, validation) per epoch — which is how the committed training figure of the documentation is produced without anything being fitted at build time.

assemble_surrogate(net, spec, box, scaling, provenance) -> NeuralSurrogate

Bundle a trained network with the box it was trained on.

Centralized on purpose: the surrogate's validity limits are the sampling box, and the feature list is the box's names. Letting a caller pass them separately is how the two drift apart.

network_widths(opts, box, spec) -> Vector{Int}

The full [n_features, hidden…, n_outputs] architecture implied by the feature box and the output specification. Both the fallback and the extension go through this, so the network the extension trains is the network the surrogate expects.

MLP(layers...)

Feed-forward stack of NNDense layers, callable on a feature vector.

The layers are held in a Tuple, so the forward pass is fully inferred and a small network costs one allocation per layer and nothing else.

NNDense(W, b, σ)

One fully connected layer, x ↦ σ.(W * x + b).

Named NNDense rather than Dense on purpose: the training extension has Lux.Dense in scope, and two Dense types in one file is a trap.

glorot_mlp(rng, widths; hidden = :tanh, output = :identity) -> MLP

Fresh network with Glorot-uniform weights and zero biases: layer k draws from $\pm\sqrt{6/(n_\mathrm{in}+n_\mathrm{out})}$, the scaling that keeps the forward variance roughly constant through a tanh stack.

widths is [n_in, h₁, …, n_out]. Every hidden layer takes the hidden activation; the output layer takes output, which should stay :identity — the physical range of a Hill-tensor component is handled by the output transform of the surrogate, not by squashing the last layer.

softplus(x)

Smooth positive activation $\log(1 + e^x)$, evaluated in the numerically stable form $\log(1 + e^{-|x|}) + \max(x, 0)$ so that large |x| neither overflows nor loses the linear branch.

activation(name::Symbol) -> Function

Look up an activation by the name used in a serialized surrogate. Raises on an unknown name rather than silently substituting a default, because a wrong activation is a silently wrong tensor.

activation_name(σ) -> Symbol

Reverse lookup, for serialization.

layer_widths(m::MLP) -> Vector{Int}

The [n_in, h₁, …, n_out] description of the architecture — what the training extension needs to build an isomorphic Lux.Chain.

layer_activations(m::MLP) -> Vector{Symbol}

Activation name of each layer, output layer included.

nparams(m::MLP) -> Int

Total number of weights and biases.

save_surrogate(path, s::NeuralSurrogate) -> String

Write s to path as pretty-printed JSON and return the path.

The round trip is exact: weights go out as Float64 decimals with full precision, so load_surrogate(save_surrogate(p, s)) reproduces every prediction bit-for-bit.

load_surrogate(path) -> NeuralSurrogate

Read a surrogate written by save_surrogate.

Needs neither Lux nor any other training dependency: a model is trained once and evaluated everywhere.

model_path(name) -> String

Absolute path of a shipped model, with or without the .json extension:

load_surrogate(model_path("spheroid_hill_iso_elastic"))

Lists what is available when the name is not found, which is more useful than a bare SystemError from the reader.

shipped_models() -> Vector{String}

Names of the surrogates committed with the package, loadable by model_path.

SURROGATE_FORMAT

Version of the on-disk surrogate format. Bumped when a field changes meaning; load_surrogate refuses a newer major version rather than guessing.

MODEL_DIR

Directory of the surrogates shipped with the package. See model_path.