API — Localization & contribution

Localization tensors

MeanFieldHomogenization.Core.strain_strain_loc Function

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

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

MeanFieldHomogenization.Core.stress_strain_loc Function

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

MeanFieldHomogenization.Core.strain_stress_loc Function

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

MeanFieldHomogenization.Core.stress_stress_loc Function

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

MeanFieldHomogenization.Core.gradient_gradient_loc Function

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

MeanFieldHomogenization.Core.flux_gradient_loc Function

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

MeanFieldHomogenization.Core.gradient_flux_loc Function

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

MeanFieldHomogenization.Core.flux_flux_loc Function

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

Contribution tensors

MeanFieldHomogenization.Core.stiffness_contribution Function

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

MeanFieldHomogenization.Core.compliance_contribution Function

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

MeanFieldHomogenization.Core.conductivity_contribution Function

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.

MeanFieldHomogenization.Core.resistivity_contribution Function

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

The amount × contribution seam

MeanFieldHomogenization.Core.delta_stiffness Function

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.

MeanFieldHomogenization.Core.delta_compliance Function

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.

MeanFieldHomogenization.Core.delta_conductivity Function

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

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

MeanFieldHomogenization.Core.delta_resistivity Function

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.

For a flat inclusion the four three-argument seams share one geometric prefactor, crack_density_factor.

Inclusion traits and bundled seams

MeanFieldHomogenization.Core.is_homogeneous_inclusion Function

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(::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.

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.

MeanFieldHomogenization.Core.loc_and_stiffness Function

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.

Contract

The returned pair must be bitwise identical to (strain_strain_loc(incl, P₁, P₀; kw...), stiffness_contribution(incl, P₁, P₀; kw...)). The generic fallback is that pair; specializations are a pure performance concern and must not reassociate the arithmetic.

Internal seam — not exported.

MeanFieldHomogenization.Core.loc_and_stress_average Function

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.

MeanFieldHomogenization.Core.compliance_and_stiffness_contribution Function

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.

Custom (user-defined) inclusions

See Custom inclusions for the tutorial and Adding a new inclusion for the full contract.

MeanFieldHomogenization.CustomInclusions Module

MeanFieldHomogenization.CustomInclusions

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

Two exports:

• CustomInclusion — a concrete, callback-driven inclusion, the value-level analog 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.

MeanFieldHomogenization.Core.AbstractCustomInclusion Type

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

MeanFieldHomogenization.CustomInclusions.CustomInclusion Type

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

Option	Default	Meaning
semi_axes	nothing	half-dimensions of an equivalent ellipsoidal envelope, if your morphology has one. Purely descriptive: no kernel reads it, and it only serves shape_tensor. Leave it out for a shape that has no such envelope.
dim	3	spatial dimension, when semi_axes is not given
basis	from euler_angles	local frame; for a flat object column 3 is the normal
euler_angles	()	ZYZ angles, ignored when basis is given
homogeneous	true	value returned by is_homogeneous_inclusion
density_factor	nothing	prefactor of the amount × contribution seam; set it (e.g. 4π/3) to register the phase with a CrackDensity amount

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	Callbacks	Signature
A — Hill tensor	hill_tensor	(P₀; kw...) -> Tens
B — localization	strain_strain_loc (+ stress_strain_loc if heterogeneous), gradient_gradient_loc (+ flux_gradient_loc)	(P₁, P₀; kw...) -> Tens
C — contribution, solid	stiffness_contribution, compliance_contribution, conductivity_contribution, resistivity_contribution	(P₁, P₀; kw...) -> Tens
C — contribution, flat	same names, with density_factor set	(P₀; kw...) -> Tens

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.

Gate B and heterogeneous inclusions

The strain-side localization determines the stress-side one only through A_σε = C₁ : A_εε, which needs a single uniform C₁. If you set homogeneous = false, supply stress_strain_loc (and flux_gradient_loc in transport) as well — otherwise the average stress in the inclusion is silently wrong, and with it SelfConsistent and AsymmetricSelfConsistent, the two schemes that consume it. Dilute and MoriTanaka need only (A, N) and stay correct, which is precisely what makes the omission easy to miss.

A heterogeneous inclusion also has no single property to feed the Voigt and Reuss bounds — they need its internal volume fractions — so those are unavailable unless the type supplies a layer-wise average (Schemes._layer_voigt / Schemes._layer_reuss, plus Schemes.has_layer_average). Every scheme that consumes localization or contribution tensors is unaffected, AsymmetricSelfConsistent included.

Gate C and the dilute concentration tensor

Contribution tensors alone do not determine A. A volume-fraction phase entered through gate C therefore works with Voigt, Reuss, Dilute, DiluteDual, Maxwell, PonteCastanedaWillis and DifferentialScheme, but not with MoriTanaka, SelfConsistent or AsymmetricSelfConsistent, whose kernels also need A. Density-based (flat) phases are unaffected: those kernels reconstruct A from ℍ:C₀.

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()
add_phase!(rve, :m, Ellipsoid(1.0), Dict(:C => iso_stiffness(10.0, 6.0)); fraction = :rest)
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.

MeanFieldHomogenization.CustomInclusions.CustomShape Type

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.

MeanFieldHomogenization.CustomInclusions.check_inclusion_interface Function

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.

Finite-element inclusions

Requires a finite-element backend: Ferrite, FerriteGmsh and Gmsh serve every morphology, while Gridap and GridapGmsh serve the crack and the axisymmetric ones. Five morphologies, one method — the same corrected finite Eshelby cell and the same entry into the schemes: Finite-element inclusions for the elliptical crack, A recycled-concrete aggregate for the sphere with an off-center core, Concave pores for the superspherical and superspheroidal cavities, and A layered spheroid, meshed for N nested spheroids.

MeanFieldHomogenization.FiniteElements Module

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

Type	Morphology	Discretization	Entry gate
FEEllipticCrack	flat elliptical crack	3-D tetrahedra	COD tensor → crack algebra
FEExcenteredSphere	sphere with an off-center spherical core	axisymmetric Fourier	B — the two localization tensors
FESupershapePore	superspherical cavity, cube-symmetric	3-D tetrahedra, whole cell or one octant	B — the strain side alone, the stress side being zero
FEAxiSupershapePore	superspheroidal cavity, a solid of revolution	axisymmetric Fourier	B — the strain side alone
FEAxiLayeredSpheroid	N nested coaxial spheroids, free semi-axes	axisymmetric Fourier	B — both localization tensors

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 generics of FEBackend, which come in three groups: the crack, the axisymmetric contract (ten methods, the tenth being fe_axi_pore_boundary, needed by a cavity alone since it has no interior to average over) and the three-dimensional cell.

Two backends exist: FerriteBackend (import Ferrite, FerriteGmsh, Gmsh), which implements all three groups, and GridapBackend (import Gridap, GridapGmsh), which implements the crack and the axisymmetric ones. 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, docs/src/applications/recycled_aggregate.md and docs/src/applications/concave_pores.md.

MeanFieldHomogenization.FiniteElements.FECache Type

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.

MeanFieldHomogenization.FiniteElements.fe_assembly_count Function

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.

MeanFieldHomogenization.FiniteElements.fe_reset! Function

fe_reset!(incl) -> incl

Drop the cached discretization and every memoized response tensor.

Choosing a backend

MeanFieldHomogenization.FiniteElements.FEBackend Type

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.

MeanFieldHomogenization.FiniteElements.AutoBackend Type

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.

MeanFieldHomogenization.FiniteElements.FerriteBackend Type

FerriteBackend()

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

MeanFieldHomogenization.FiniteElements.GridapBackend Type

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.

Writing a backend

A backend is sixteen methods and nothing else — nine for the axisymmetric solve, seven for the crack. The Fourier operators, the boundary data, the fixed point of the corrected boundary condition and the memoization are shared, and the driver closes the strain operator and the azimuthal projection over the mode before handing them over, so an implementation never sees a Fourier mode or a physics: only "this many scalar fields, this operator, this projection".

ext/MeanFieldHomogenizationGridapExt/ is the shorter of the two implementations and the one to read first.

The axisymmetric solve

MeanFieldHomogenization.FiniteElements._build_gmsh_axi_model Function

_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 centered 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().

MeanFieldHomogenization.FiniteElements.fe_axi_grid Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_grid_counts Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_region_volume Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_mode Function

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

MeanFieldHomogenization.FiniteElements.fe_axi_dof_split Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_set_dirichlet! Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_stiffness Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_average Function

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.

MeanFieldHomogenization.FiniteElements._resolve_backend Function

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

The crack

MeanFieldHomogenization.FiniteElements._build_gmsh_crack_model Function

_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 centered 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().

MeanFieldHomogenization.FiniteElements._weld_msh_crack_front Function

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

MeanFieldHomogenization.FiniteElements.fe_crack_grid Function

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.

MeanFieldHomogenization.FiniteElements.fe_crack_counts Function

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.

MeanFieldHomogenization.FiniteElements.fe_crack_space Function

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.

MeanFieldHomogenization.FiniteElements.fe_crack_dof_split Function

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.

MeanFieldHomogenization.FiniteElements.fe_crack_set_dirichlet! Function

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.

MeanFieldHomogenization.FiniteElements.fe_crack_stiffness Function

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.

MeanFieldHomogenization.FiniteElements.fe_crack_mean_jump Function

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.

The three-dimensional cell

Eight methods for the third cell family. Two things distinguish it from the crack's: the geometry is quadratic, the mid-edge nodes of the inclusion boundary having been moved onto the exact shape, and one contract serves two physics — a scalar temperature and a vector displacement — with only the material and the two averages differing.

MeanFieldHomogenization.FiniteElements._build_gmsh_cell_model Function

_build_gmsh_cell_model(gmsh, shape, opts) -> NamedTuple

Build and mesh the cell in the already initialized gmsh session. Returns the two surfaces actually used, the radius, the applied sizes and the snapping report.

The meshed domain is the matrix shell between the inclusion surface and the outer sphere: a pore is not meshed at all, its boundary being naturally flux-free in conduction and traction-free in elasticity. That is why the volume is built from two surface loops with the outer one first — the inclusion is a hole.

cavity_volume is the volume of the body the solve actually sees, measured on the curved boundary once it has been snapped. Normalizing a localization tensor by the exact volume while solving on the meshed one would fold a geometry error into the answer with nothing to reveal it; reporting both is what makes that error visible instead.

MeanFieldHomogenization.FiniteElements._snap_cell_surface_to_shape! Function

_snap_cell_surface_to_shape!(gmsh, shape, tag; tol, passes) -> NamedTuple

Move the mid-edge nodes of a second-order surface onto the exact shape, backing off wherever that would invert a neighboring element.

gmsh.model.mesh.setOrder(2) places every mid-edge node at the midpoint of the straight segment. It cannot do better: the surface was handed to it as a discrete entity, so it has no idea what surface the node is supposed to be on. Projecting those nodes radially — exact and single-valued here, the body being star-shaped — turns the boundary into a genuine second-order approximation and takes the geometry error from to . Measured on this family it is worth about a factor of 16 in cost at equal geometric accuracy.

Why it is safeguarded

The displacement is , so asymptotically nothing can be inverted. That argument fails on a concave shape, where the curvature near the conical points on the axes is not small against — and, unlike an ordinary under-resolution, it never becomes small, because at a conical point the curvature is unbounded.

So each node carries a blend factor. All are snapped fully, the element qualities are read back, and any node belonging to an element below tol has its factor halved; repeat. limited counts the nodes that had to back off and min_blend says how far the worst one did.

Refinement does not remove it, and it is not meant to. Measured on  , radius_ratio = 3:

level	tetrahedra	limited	fraction	min_blend
2	5 604	6 / 258	2.3 %	0.5
3	21 736	37 / 1026	3.6 %	0.25
4	102 641	73 / 4098	1.8 %	0.125

The count grows and the worst blend halves at every level, because the same non-differentiable feature is being approached by ever-smaller elements whose curvature relative to it stays infinite. What is stable is the fraction, a few percent, and that is the number to watch: a large fraction means the level really is too coarse for that , a few percent means the conical points are being handled as they must be.

What separates a shape that needs the safeguard from one that does not is unbounded curvature, not sharpness. At level 3 the sphere, the mildly concave  , and the convex   all report limited = 0; only   limits. And the octahedron,  , reports maxd = 1.2e-16: it has edges and vertices everywhere, yet there is nothing to snap at all, because its faces are flat and the mid-edge nodes are already exactly on the surface. Edges are harmless; conical points are not.

The nodes are shared with the adjacent tetrahedra, so setNode moves the volume mesh along with the surface. That is the point, not a side effect.

MeanFieldHomogenization.FiniteElements._snap_cell_surface_to_sphere! Function

_snap_cell_surface_to_sphere!(gmsh, R, tag) -> NamedTuple

The counterpart of _snap_cell_surface_to_shape! for the outer boundary: push every node of a second-order surface onto the sphere of radius R. No safeguard is needed — the outer mesh is uniform and far from anything sharp.

MeanFieldHomogenization.FiniteElements.fe_cell_grid Function

fe_cell_grid(backend, shape, opts)

Backend-native mesh of the cell around shape, carrying the cell set "matrix" and the boundary sets "inclusion" and "outer" of _build_gmsh_cell_model.

The mesh is second order with a curved boundary: the mid-edge nodes of the inclusion surface have been moved onto the exact shape, so the geometric interpolation must be quadratic too. A backend that silently used a linear geometry here would throw away the whole gain and report a plausible wrong answer.

MeanFieldHomogenization.FiniteElements.fe_cell_counts Function

fe_cell_counts(backend, grid) -> (; ncells, nnodes, area_inclusion, area_outer)

Mesh diagnostics: cell and node counts, and the area of each boundary set. The outer area is compared against — it is what says the outer surface was snapped onto the sphere rather than left a polyhedron.

MeanFieldHomogenization.FiniteElements.fe_cell_space Function

fe_cell_space(backend, grid, order, ncomp) -> space

ncomp scalar Lagrange fields of degree order on grid — ncomp = 1 for conduction, 3 for elasticity — with a quadrature exact to degree 2 * order, a facet quadrature of the same degree, and a quadratic geometric interpolation to follow the curved boundary.

The dof numbering must span the whole space, with no Dirichlet elimination: the driver splits the dofs itself so that one factorization serves all twelve right-hand sides.

MeanFieldHomogenization.FiniteElements.fe_cell_dof_split Function

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

Total dof count and the two index vectors, presc being the dofs of the "outer" boundary and nothing else. The inclusion surface carries no condition at all — that is what makes it a pore.

MeanFieldHomogenization.FiniteElements.fe_cell_set_dirichlet! Function

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

Write u[d] = f(x_d)[k] for every dof d of the "outer" boundary, x_d being its node and k its component; for a scalar field f returns a number.

Called once per right-hand side and must not touch the matrix: the driver assembles and factorizes once, then calls this twelve times.

MeanFieldHomogenization.FiniteElements.fe_cell_stiffness Function

fe_cell_stiffness(backend, space, material) -> AbstractMatrix

Stiffness over the whole dof numbering, over the matrix region:

conduction   K = ∫ ∇v ⋅ 𝐊 ⋅ ∇u dΩ         `material` a 3×3 matrix
elasticity   K = ∫ ε(v) : ℂ : ε(u) dΩ      `material` a SymmetricTensor{4,3}

Two methods on one name, told apart by the type of material. Both are isotropic: the corrected boundary condition uses a closed-form dipole field, and the driver refuses anything else long before reaching here, so a backend may work from the scalar or from (λ, μ) instead of from the full object.

MeanFieldHomogenization.FiniteElements.fe_cell_mean_gradient Function

fe_cell_mean_gradient(backend, space, u, V) -> NTuple{3}

over the inclusion, as the surface integral   , with V the volume of the inclusion.

The inclusion is not meshed, so its average is reached through the divergence theorem on its boundary — which is also why it costs nothing beyond the solve.

On the sign. The normal a backend reports on a facet set points out of the meshed region, hence into the pore, so it is the opposite of the inclusion's own outward normal. The minus sign belongs here, and it is the same one the crack and axisymmetric drivers carry in two different places for the same reason. Getting it wrong is not invisible: the corrected answer then carries exactly twice the truncation bias instead of none.

MeanFieldHomogenization.FiniteElements.fe_cell_mean_strain Function

fe_cell_mean_strain(backend, space, u, V) -> NTuple{6}

over the inclusion in Kelvin-Mandel, as   .

Same normal convention, and the same minus sign, as fe_cell_mean_gradient.

MeanFieldHomogenization.FiniteElements._cell_close_dipole Function

_cell_close_dipole(L_s, L_u, F) -> (L, dipole_norm)

Superpose the remote and unit-dipole families: L = (𝕀 − L_u F)⁻¹ L_s.

dipole_norm = ‖L_u F‖ is the diagnostic to watch rather than the answer: the dipole term is , so a log-log slope of against R says the correction is doing what it claims.

MeanFieldHomogenization.FiniteElements._cell_conduction_localization Function

_cell_conduction_localization(backend, space, k₀, V; sign = -1)
    -> (; A, A_uncorrected, L_s, L_u, dipole_norm)

The gradient localization tensor of the cavity, corrected and uncorrected, as a 3×3 in the cell's frame.

Six solves on one factorization: three with the remote gradient   and three with a unit-moment dipole. sign selects the sign of the polarization factor and exists so that a test can flip it and watch the factor of two appear.

MeanFieldHomogenization.FiniteElements._cell_elastic_localization Function

_cell_elastic_localization(backend, space, C₀, μ, ν, V; sign = -1)
    -> (; A, A_uncorrected, L_s, L_u, dipole_norm)

The strain localization tensor of the cavity, corrected and uncorrected, as a 6×6 in Kelvin-Mandel.

Twelve solves on one factorization: six with the remote strain   of the Kelvin basis — which is what makes a unit load a unit load, the map being an isometry — and six with a unit-moment dipole.

MeanFieldHomogenization.FiniteElements._cell_kelvin_basis Function

_cell_kelvin_basis(i) -> 3×3

The i-th Kelvin-Mandel basis tensor, orthonormal for the Frobenius inner product — which is what makes a unit load in this basis a unit load, and lets the same six matrices serve both as remote strains and as dipole moments.

MeanFieldHomogenization.FiniteElements._cell_outer_radius Function

_cell_outer_radius(shape, opts) -> R

radius_ratio times the shape's bounding_radius — see FECellMeshOptions for why it is not a.

MeanFieldHomogenization.FiniteElements._cell_solver Function

_cell_solver(backend, space, K)

Close over one factorization of the free-free block so that every right-hand side is a back-substitution. The prescribed dofs are eliminated by hand rather than applied into the matrix, which would consume it.

Superspherical and superspheroidal pore

The inclusion type. It enters as a heterogeneous inclusion, which for a cavity is the truthful answer rather than a convenience: is_homogeneous_inclusion asks whether a single describes the interior, and inv(0) is meaningless. The package's exact identities then take over, and with a cavity's stress-side localization being identically zero they collapse to     and    — the right answer for a pore, with no anywhere in it. So the phase property handed to add_phase! is genuinely ignored, and no contribution tensor is overridden.

MeanFieldHomogenization.FiniteElements.FESupershapePore Type

FESupershapePore(shape; opts = FECellMeshOptions(), basis = nothing,
                 euler_angles = (), backend = AutoBackend())

A cavity of superspherical or superspheroidal shape, whose response is computed on the truncated cell with a corrected boundary condition.

shape is a Supersphere or a Superspheroid — see those for the convention. The phase property handed to add_phase! is a placeholder and is ignored: the solve is the exact cavity problem, and there is no in the answer. Pass the matrix's own stiffness and nothing is lost.

Both physics are served by the same object and the same mesh — a scalar solve for conduction, a vector one for elasticity — so one type covers homogenize(rve, scheme, :C) and homogenize(rve, scheme, :K).

The reference medium must be isotropic: the corrected boundary condition uses the closed-form Kelvin dipole field. Under an iterative scheme, whose running estimate for a cubic inclusion is cubic and not isotropic, add symmetrize = IsoSymmetrize() to the phase — the error message says so rather than letting the scheme converge to something wrong.

Requires a finite-element backend: import Ferrite, FerriteGmsh, Gmsh.

Symmetry

A supersphere in an isotropic matrix is cubic, so its compliance contribution has three independent constants and not two. A superspheroid is transversely isotropic about . Neither class is projected onto by default: the raw tensors come out, and cubic_residual is the free error bar on the cubic case — the distance to a class the answer belongs to by group theory is discretization error and nothing else.

Example

using MeanFieldHomogenization
import Ferrite, FerriteGmsh, Gmsh

pore = FESupershapePore(Supersphere(1.0, 0.6);
                        opts = FECellMeshOptions(; level = 3, radius_ratio = 3.0))

C₀ = iso_stiffness(0.8333, 0.3846)          # E = 1, ν = 0.3
rve = RVE()
add_phase!(rve, :m, Ellipsoid(1.0), Dict(:C => C₀); fraction = :rest)
add_phase!(rve, :pores, pore, Dict(:C => C₀); fraction = 0.05)
homogenize(rve, MoriTanaka(), :C)

See also FECellMeshOptions, fe_cell_localization, fe_cell_mesh_report.

MeanFieldHomogenization.FiniteElements.SupershapePoreShape Type

SupershapePoreShape

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

MeanFieldHomogenization.FiniteElements.fe_cell_localization Function

fe_cell_localization(pore, P₀) -> A

The localization tensor of the cavity — for a fourth-order reference, for a second-order one — in the global frame.

Memoized on the reference medium expressed in the pore's own frame, which has a useful consequence: a whole family of orientations of the same shape in the same matrix shares one finite-element resolution. So a one-shot scheme costs a single solve even under an orientation average.

MeanFieldHomogenization.FiniteElements.fe_cell_mesh_report Function

fe_cell_mesh_report(pore) -> NamedTuple

What the mesh actually is: cell and node counts, the two boundary areas, the radius, the applied sizes, the snapping report, and the three volumes worth comparing — the closed form, the flat triangulation it was built from, and the curved boundary the solve sees.

Those three are the geometry error, made visible. The gap between the flat and the curved volume is what snapping the mid-edge nodes bought; the gap between the curved volume and the closed form is what is left.

MeanFieldHomogenization.FiniteElements.has_surrogate Function

has_surrogate(pore, order) -> Bool

Whether this pore answers order-order physics from a network rather than from a mesh — 4 elasticity, 2 transport.

It is the question that decides whether the response can be differentiated with respect to the morphology. A finite-element solve cannot be: it runs in Float64 and memoizes on the reference medium alone, so the derivative would come back a silent zero and the package refuses it outright. A surrogate can, and for a shape family indexed by a single parameter p that is most of the reason to train one.

MeanFieldHomogenization.FiniteElements.pore_shape_params Function

pore_shape_params(pore) -> NamedTuple

The morphology parameters by name: (; a, p) for a supersphere, (; a, c, p) for a superspheroid.

One list serving two purposes, which is why it is a function rather than a convention written twice: it is what a surrogate names as its features, and it is what the sensitivity API differentiates with respect to.

MeanFieldHomogenization.FiniteElements._pore_surrogate_response Function

_pore_surrogate_response(surrogate, pore, P₀) -> AbstractTens

Evaluate a trained surrogate in place of the finite-element solve. Declared here and implemented in NeuralInclusions, the only module that knows what a NeuralSurrogate is and which is included later.

MeanFieldHomogenization.FiniteElements._rebuild_pore_shape Function

_rebuild_pore_shape(shape, Val(name), value) -> shape

The shape with one named parameter replaced. What makes a surrogate-backed pore differentiable in its own morphology: the sensitivity API rebuilds the geometry with a ForwardDiff.Dual in place of one parameter, and everything downstream promotes.

MeanFieldHomogenization.Schemes._geom_field Function

_geom_field(geom, name) -> value

Read a named geometry parameter. getfield by default, and overridable, which is what lets a geometry expose parameters that are not literal fields — a type holding a shape object rather than its scalars, say. The reader and the rebuilder, _replace_geom_field, must agree on the same name set or a derivative silently differentiates the wrong thing.

MeanFieldHomogenization.FiniteElements.FEAxiSupershapePore Type

FEAxiSupershapePore(shape; opts, basis, euler_angles, backend,
                    elastic = nothing, transport = nothing, guard = :warn)

A superspheroidal cavity solved by axisymmetric Fourier finite elements.

shape is a Superspheroid; its axis of revolution is the third axis of basis. The answer is transversely isotropic — five constants in elasticity, two in transport — which is exactly what the Fourier modes deliver: mode 0 gives a 2×2 block, modes 1 and 2 a scalar each.

elastic and transport take a trained surrogate in place of the solve, one per physics, exactly as for FESupershapePore. That is the only route to a sensitivity with respect to the morphology: the finite-element solve runs in Float64 and memoizes on the reference medium alone, so a derivative in p would come back a silent zero.

radius_ratio multiplies the shape's bounding radius, not a — for an elongated superspheroid the two differ by the aspect ratio, and using a lets the outer boundary come within 1.2c of the body.

MeanFieldHomogenization.FiniteElements.FEAxiLayeredSpheroid Type

FEAxiLayeredSpheroid(axis_radii, disk_radii, props; kwargs...)

An N-layer coaxial, concentric spheroid in an isotropic matrix, solved by Fourier axisymmetric finite elements on the meridian half-plane.

axis_radii[ℓ] and disk_radii[ℓ] are layer ℓ's semi-axes along and across the revolution axis, ascending, core first — the argument order of LayeredSpheroid, so confocal_layer_radii feeds both. props is one isotropic modulus tensor per layer: Tens{4,3} for elasticity, Tens{2,3} for transport, and the order chosen fixes the object's physics.

keyword	default	meaning
interfaces	all perfect	per-layer interface; only PerfectInterface is implemented
opts	FEAxiMeshOptions()	mesh density and cell radius
basis	canonical	the revolution axis is the basis's third column
backend	AutoBackend()	finite-element backend

Enters through gate B, with both localization tensors measured on the same solve: the inclusion is heterogeneous, so the stress side is not ℂ₁ : 𝔸_εε for any single ℂ₁.

Examples

# The validation slice: confocal, described exactly as the analytic type is.
ar, dr = confocal_layer_radii(2.0, 1.0, (0.3, 0.7))
K = (TensISO{3}(5.0), TensISO{3}(2.0))
fe  = FEAxiLayeredSpheroid(ar, dr, K)
ana = LayeredSpheroid(ar, dr, K)

# Free radii: an oblate core in a prolate shell. No closed form covers it.
FEAxiLayeredSpheroid((0.4, 1.4), (0.9, 1.0), K)

MeanFieldHomogenization.FiniteElements.LayeredSpheroidShape Type

LayeredSpheroidShape

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

MeanFieldHomogenization.FiniteElements.check_nested_spheroids Function

check_nested_spheroids(axis_radii, disk_radii, R = Inf)

Verify that the layers are strictly nested and inside the cell, and throw a message naming the offending layer otherwise.

Two coaxial concentric ellipses are nested if and only if both semi-axes of the inner one are no larger than the outer's, so the whole condition is two comparisons per layer. It is worth checking rather than trusting: a violation does not produce a wrong answer, it produces a self-intersecting geometry that gmsh rejects from deep inside its own meshing pipeline, with a message that names neither the layer nor the semi-axis.

disk_radii are the transverse semi-axes and axis_radii those along the revolution axis — the argument order LayeredSpheroid uses.

Consecutive layers must be strictly ascending in both semi-axes. Two coinciding boundaries pass every containment argument and carry no volume, but they cannot be meshed: the element size is capped at a fraction of a layer's thickness, so a zero-thickness layer asks gmsh for a zero-sized element. A sphere is still a legitimate layer — that is a == c within one layer, which is unrestricted, and it is the slice where LayeredSphere is the closed form.

MeanFieldHomogenization.FiniteElements.axi_layer_set Function

axi_layer_set(ℓ) -> String

Physical-group name of layer ℓ, counted from the core. The matrix keeps AXI_SET_MATRIX, so a driver reads N + 1 regions from one list.

MeanFieldHomogenization.FiniteElements.layer_volumes Function

layer_volumes(incl) -> Vector

Volume of each layer, 4π/3 · a² c of its own boundary minus the one inside it. Closed form for any semi-axes, confocal or not — which is what makes it the mesh check that is always available.

MeanFieldHomogenization.FiniteElements.layer_fractions Function

layer_fractions(incl) -> Vector

Each layer's share of the inclusion's volume, summing to one.

MeanFieldHomogenization.FiniteElements.AxiSupershapePoreShape Type

AxiSupershapePoreShape

Shape trait of FEAxiSupershapePore. No kernel dispatches on it — the inclusion supplies its own localization tensors.

MeanFieldHomogenization.FiniteElements.fe_axi_pore_localization Function

fe_axi_pore_localization(pore, P₀) -> Tens

The cavity's localization tensor, memoized on the reference medium expressed in the pore's own frame — so a whole family of orientations of the same shape in the same matrix shares one solve.

MeanFieldHomogenization.FiniteElements.fe_axi_pore_breakdown Function

Bypass the cache and return every intermediate, including A_uncorrected.

MeanFieldHomogenization.FiniteElements.fe_axi_pore_mesh_report Function

fe_axi_pore_mesh_report(pore) -> NamedTuple

Cell and node counts, and the one comparison that says whether the meridian spline resolves the profile: the measured revolution volume of the matrix against V_Ω - V_D in closed form. The gap is the geometry error, and it is the number to watch when p goes concave.

MeanFieldHomogenization.FiniteElements.fe_axi_pore_boundary Function

fe_axi_pore_boundary(backend, mode, u, dofmap, proj, set) -> Vector

Line integral of (ū ⊗ n)ˢ on the meridian trace of a cavity wall, projected onto the Fourier mode, with the measure ρ dl and unnormalized. n is the facet normal, outward from the matrix and therefore into the cavity, so the caller supplies the sign.

The tenth method of the axisymmetric contract, and the only one a cavity needs that a solid inclusion does not. A cavity has no interior to integrate over, so ⟨ε⟩_D has to come from its boundary — exactly as fe_cell_mean_strain does in three dimensions.

Not by the divergence identity on the matrix. That route is algebraically exact and numerically hopeless: it computes V_D as the difference of V_Ω and V_M, whose ratio is (R/a)³, so it loses two digits at R/a = 4 and more as the cell grows. Measured: the implied cavity volume is 2.9 % off at V_Ω/V_D = 8 and 9.2 % off at 216, and the localization error tracks it.

dofmap is _axi_dof_map(m), needed here and not in fe_axi_stiffness because Bop already folds it in: the boundary term uses the raw (ū_ρ, ū_θ, ū_z) while the solved unknowns are the mapped ones.

The amplitude is returned in the same Kelvin ordering as _axi_B_elast produces — (ρρ, θθ, zz, √2 θz, √2 ρz, √2 ρθ) — so the same proj applies. Its θθ entry is identically zero, the facet normal having no azimuthal component; the metric term that ε_θθ carries has no counterpart in a tensor product, and the trace identity u_ρn_ρ + u_θn_θ = u₁n₁ + u₂n₂ is what makes that consistent.

MeanFieldHomogenization.FiniteElements._superspheroid_meridian Function

_superspheroid_meridian(shape, n) -> Vector{NTuple{2, Float64}}

n points of the meridian profile from the equator (a, 0) to the pole (0, c), both included, as (ρ, z) pairs.

Exact by construction: the parametrization ρ = a t^{1/m}, z = c(1−t)^{1/m} satisfies (ρ/a)^m + (z/c)^m = t + (1−t) = 1 for every t, so the points lie on the surface to round-off and level_set returns zero on each — which the tests check rather than assume.

The step clusters at both ends, through t = (1 − cos πs)/2. Uniform in t would put its points where the profile is flat and starve the two corners, whose curvature is unbounded for a concave shape.

MeanFieldHomogenization.FiniteElements._fe_frame Function

_fe_frame(incl) -> 3×3

Rotation matrix whose columns are the inclusion's local axes in global coordinates. Shared by every finite-element inclusion: each one solves in its own frame and rotates the result out.

Elliptical crack (3-D)

MeanFieldHomogenization.FiniteElements.FEEllipticCrack Type

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 MeanFieldHomogenization, 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()
add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => C₀); fraction = :rest)
add_phase!(rve, :cracks, crack, Dict(:C => C₀); density = 0.05,
           symmetrize = IsoSymmetrize())
homogenize(rve, MoriTanaka(), :C)

See also FEMeshOptions, FECache.

MeanFieldHomogenization.FiniteElements.FEMeshOptions Type

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

Discretization settings of a finite-element inclusion.

Field	Default	Meaning
radius_ratio	5.0	radius of the surrounding ball of matrix, in units of the largest semi-axis. 5 is enough because the boundary condition is corrected — the uncorrected problem needs 10 to 40.
htipdiv	12.0	element size at the crack front, as b / htipdiv. The mesh coarsens to R/3 at the outer boundary.
order	2	polynomial order of the displacement interpolation (1 or 2). Tetrahedra are geometrically straight; order = 2 is therefore subparametric.

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

MeanFieldHomogenization.FiniteElements.fe_mesh_report Function

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 equaling πab is what says the Crack plugin split the surface cleanly and the front weld did not glue the lips back together.

MeanFieldHomogenization.FiniteElements.fe_cod_breakdown Function

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.

Sphere with an off-center core (axisymmetric Fourier)

MeanFieldHomogenization.FiniteElements.FEExcenteredSphere Type

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 center, resolved by axisymmetric Fourier finite elements.

Geometry

Symbol	Definition
a	radius of the whole inclusion
w = core_fraction	volume fraction of the core within the inclusion
a_core = a·w^(1/3)	core radius, fixed by w
α = eccentricity	offset of the core center, as a fraction of the largest offset that keeps the core inside: d = α·(a − a_core)

α = 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 MeanFieldHomogenization
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()
add_phase!(rve, :m, Ellipsoid(1.0), Dict(:C => C₀); fraction = :rest)
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.

MeanFieldHomogenization.FiniteElements.FEAxiMeshOptions Type

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

Discretization settings of an axisymmetric finite-element inclusion.

Field	Default	Meaning
radius_ratio	4.0	radius R of the surrounding ball of matrix, in units of the inclusion radius. 4 suffices because the boundary condition is corrected — the uncorrected problem needs 10 or more (Adessina et al. 2017, Table 2).
nradial	24	element size inside the inclusion, as a / nradial.
coarsening	6.0	ratio of the element size at the outer boundary to the size inside the inclusion.
order	2	polynomial order of the displacement / temperature interpolation (1 or 2).
nprofile	61	points sampled on a quarter of the meridian profile — FEAxiSupershapePore only.
tip_refine	16	element size at the profile's two corners, as a divisor of h_in — the same pore only.

Grading and nradial compound

On a concave profile the two knobs multiply rather than add: at p = 0.20 and R/a = 6, nradial = 20 meshes 68 000 cells, nradial = 28 130 000 and nradial = 40 267 000 — the last exceeding a 6 GB budget in elasticity. Refine radius_ratio before nradial: what remains after the dipole correction is truncation, and nradial = 20 is already past its plateau (going to 28 leaves H₁₁₁₁ and R₁₁ unchanged to six figures).

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.

MeanFieldHomogenization.FiniteElements.fe_axi_localization Function

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.

MeanFieldHomogenization.FiniteElements.fe_axi_breakdown Function

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.

Bypass the cache and return every intermediate, including the uncorrected pair.

MeanFieldHomogenization.FiniteElements.fe_axi_mesh_report Function

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.

fe_axi_mesh_report(incl::FEAxiLayeredSpheroid) -> NamedTuple

Mesh diagnostics: cell and node counts, and each layer's measured volume of revolution against its closed form.

The closed form is 4πa²c/3 of the layer's own boundary minus the one inside it, which exists for any semi-axes — confocal or not. So this is the check that remains available on a geometry with no analytic solution at all, and it is worth running before a long solve on a nest that has not been tried: the volume_error column catches a mesh that failed to resolve a thin layer, which otherwise looks perfectly reasonable.

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.

MeanFieldHomogenization.FiniteElements.core_radius Function

core_radius(incl) -> Real

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

MeanFieldHomogenization.FiniteElements.core_offset Function

core_offset(incl) -> Real

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

MeanFieldHomogenization.FiniteElements.tensor_order Function

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.

MeanFieldHomogenization.FiniteElements.ExcenteredSphereShape Type

ExcenteredSphereShape

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

The three-dimensional cell around a non-ellipsoidal shape

The third cell family. The inclusion surface has no CAD representation, so it is built analytically by shape_surface and handed to gmsh as a discrete entity; the interior size follows an exact radial law, the cell being star-shaped about its center.

Two volume measures rather than one, and the difference is the point: the flat mesh_volume cannot see a curved boundary at all, so fe_cell_curved_volume is what shows what snapping the mid-edge nodes onto the exact shape actually bought.

MeanFieldHomogenization.FiniteElements.FECellMeshOptions Type

FECellMeshOptions(; kwargs...)

Geometry and resolution of the truncated cell around a supershape. Every size is a knob: the machine budget is a hard constraint, and the useful resolution is not known before a convergence study.

Keyword	Default	Meaning
radius_ratio	4.0	outer radius as a multiple of the shape's bounding radius
level	4	subdivision level of the inclusion surface:   triangles
outer_level	3	subdivision level of the outer sphere
relax	60	tangential relaxation sweeps on the inclusion surface
h_in, h_out	nothing	interior element size at each boundary; derived from the surface meshes when nothing
grading	1.5	exponent of the radial size law
order	2	element order; 2 is required for a curved boundary
algorithm3d	1	gmsh Mesh.Algorithm3D: 1 Delaunay (robust with a fixed boundary), 10 HXT (faster)
optimize	true	run gmsh's tetrahedron optimizer
max_dofs	200_000	refuse to factorize above this
min_free_gb	6.0	refuse to start a factorization with less than this available
verbose	false	let gmsh talk

Three of these deserve more than a line.

outer_level is not a free knob. At outer_level = 2 the outer boundary carries 128 curved triangles, which does not resolve the dipole term of the corrected boundary condition: the exact spherical-pore gate then stalls at  . Level 3 takes it to  .

radius_ratio multiplies the bounding radius, not a. For an elongated superspheroid the two differ by the aspect ratio, and using a lets the outer boundary come within of the body — where the exact spheroid gate falls from to  . For a supersphere with   the bounding radius is , so the literature's convention is preserved exactly where the literature uses it.

octant = true meshes one eighth of the cell, which is exact whenever the three coordinate planes are mirror planes of the shape — see has_coordinate_mirrors, a weaker condition than cubic symmetry that both shipped families satisfy. It is worth a factor of eight in degrees of freedom and more in factorization cost. The default is false because the mesh changes, so every reference number in the documentation and the regression tests would move with it; and because max_dofs then counts the dofs of this solve, the peak being one factorization out of the four (three in transport) the parity classes require.

max_dofs and min_free_gb are not performance knobs. See fe_available_gb: they exist because an out-of-memory kill takes the session with it.

MeanFieldHomogenization.FiniteElements.fe_cell_size_estimate Function

fe_cell_size_estimate(shape, opts) -> NamedTuple

Predict (; ntets, nnodes_p1, nnodes_p2, dofs_scalar, dofs_vector, …) without meshing, by integrating the radial size law:

The point is to find out that a set of options will exhaust the machine before gmsh spends ten minutes proving it. It is an estimate, not a measurement, and it over-predicts by about a factor of three: at  , level 3, radius_ratio = 3, it announces 184 000 vector degrees of freedom where the mesh actually carries 61 400. That direction is the desired failure mode for something whose job is to refuse a solve — but the factor is stated rather than called "of order one", so that nobody sizes a study from it.

MeanFieldHomogenization.FiniteElements.fe_cell_curved_volume Function

fe_cell_curved_volume(gmsh, tag; sub = 8) -> V

Volume enclosed by a second-order surface,   , evaluated by splitting every 6-node triangle into flat sub-triangles through the quadratic map and summing the signed cones from the origin.

This is what measures the gain from snapping the boundary: the flat mesh_volume cannot see a curved boundary at all, and would report the same number before and after.

MeanFieldHomogenization.FiniteElements.fe_cell_meshed_volume Function

fe_cell_meshed_volume(gmsh, tag = CELL_TAG_MATRIX) -> V

Volume of the meshed region, summed over its tetrahedra from their four corner nodes.

gmsh.model.occ.getMass is not an option: the volume lives in the geo kernel and is bounded by discrete surfaces, which OCC knows nothing about. Straight tetrahedra also mean this slightly under-reports a second-order mesh with a curved boundary, by the same that snapping buys back on the surface — so use it as a check on the topology, not as a precision measurement.

MeanFieldHomogenization.FiniteElements.fe_available_gb Function

fe_available_gb() -> Float64

MemAvailable from /proc/meminfo, in GiB, or Inf where that file is not readable — every platform other than Linux, and some containers.

Inf is the deliberate answer for "unknown": a guard that refuses to run because it cannot measure is worse than no guard.

Green function of the corrected boundary condition

MeanFieldHomogenization.Core.green_gradient_iso Function

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

Gradient of the Kelvin Green function of an isotropic elastic matrix C₀, evaluated at x ≠ 0.

With  ,   and    the Kelvin solution reads

so that

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.

MeanFieldHomogenization.Core.dipole_displacement_iso Function

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

Π has the dimension of a stress times a volume: for an inclusion of volume carrying a uniform polarization (i.e.     inside it),   .

For a symmetric Π the closed form collapses to

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.

MeanFieldHomogenization.Core.green_gradient_iso2 Function

green_gradient_iso2(K₀::TensISO{2,3}, x) -> SVector{3}
green_gradient_iso2(k₀::Number, x) -> SVector{3}

Gradient of the Green function of an isotropic conducting matrix of conductivity , evaluated at  .

With  ,

Throws a DomainError at the origin. Type-generic.

See also dipole_temperature_iso, and green_gradient_iso for the elastic case this mirrors.

MeanFieldHomogenization.Core.dipole_temperature_iso Function

dipole_temperature_iso(K₀, x, M) -> Number

Temperature at x radiated in an infinite isotropic conducting matrix by a point polarization of vector moment M,

M has the dimension of a flux times a volume: an inclusion of volume carrying a uniform polarization (that is,      inside it) radiates with   .

On the sign. It is the one that follows from the two definitions above, and it is worth measuring rather than trusting: a wrong sign leaves the corrected answer carrying exactly twice the truncation bias instead of none, which looks like a mesh that will not converge rather than like an algebra mistake. The factor of two is the signature, and it is what a test should flip on purpose to see.

See also dipole_displacement_iso.

Neural-surrogate inclusions

The fourth route into the contract: the response comes out of a trained network. See Neural-surrogate inclusions for the tutorial. Evaluating needs nothing beyond the package; training needs import Lux, Optimisers, Zygote.

MeanFieldHomogenization.NeuralInclusions Module

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

Type	Entry gate	For
NeuralHillInclusion	A — the Hill tensor	a morphology with a Hill tensor: contrast dependence and the ℂ₁ = ℂ₀ ⟹ 𝔸 = 𝕀 limit stay exact
NeuralLocalizationInclusion	B — both localization tensors	an internally heterogeneous morphology, which has no Hill tensor

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 labeling by a teacher, the one seam that changes between morphologies.

• training.jl — TrainingOptions and the fallback train_surrogate; the optimizer lives in MeanFieldHomogenizationLuxExt.

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

MeanFieldHomogenization.NeuralInclusions.NeuralHillInclusion Type

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.

Option	Meaning
semi_axes	the geometry, as a tuple; also the differentiable field
elastic	order-4 surrogate, for hill_tensor(incl, C₀)
transport	order-2 surrogate, for hill_tensor(incl, K₀)
basis / euler_angles	the local frame, as everywhere else in the package
guard	:warn (default), :error or :none — what to do when a feature falls outside the box the surrogate was trained on

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()
add_phase!(rve, :m, Ellipsoid(1.0), Dict(:C => iso_stiffness(20.0, 12.0)); fraction = :rest)
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.

MeanFieldHomogenization.NeuralInclusions.NeuralLocalizationInclusion Type

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.

Option	Meaning
strain / stress	the order-4 pair, 𝔸_εε and 𝔸_σε
gradient / flux	the order-2 pair, 𝔸_∇∇ and 𝔸_q∇
shape_params	named morphology parameters, e.g. (; eccentricity = 0.4, core_fraction = 0.5). They are the surrogate's features and the fields the sensitivity API differentiates, so every value must be a Number
fractions / properties	internal volume fractions and constituent properties; supplying them unlocks the Voigt and Reuss bounds, which a heterogeneous inclusion cannot otherwise serve. properties also backs the :log_mu_ratio_k contrast features

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.

Not exercised by the shipped models

The pilot trains gate-A surrogates, where the contrast dependence is exact. This type is the seam for the heterogeneous morphologies, and it is tested against an equivalent gate-A inclusion; no trained model ships for it yet.

MeanFieldHomogenization.NeuralInclusions.NeuralShape Type

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.

MeanFieldHomogenization.NeuralInclusions.StrainLocTI Type

Transversely isotropic 4th-order tensor without major symmetry: the six Walpole components (ℓ₁, …, ℓ₆) of TensTI{4,·,6}, for the strain-side localization tensor 𝔸_εε of a morphology axisymmetric about n.

Six and not five, because a localization tensor is genuinely not major-symmetric: 𝔸_εε = [𝕀 + ℙ:(ℂ₁−ℂ₀)]⁻¹ is the inverse of a product of two major-symmetric tensors, which does not commute. For an oblate spheroid at a contrast of 2 the defect reaches 10 % and ℓ₃ ≠ ℓ₄ outright, so projecting onto the five-component form would lose a few percent.

Six and not eight, because TensTI{4,·,8} also carries the couplings ℓ₇, ℓ₈ that are antisymmetric in an index pair. A tensor mapping symmetric strains to symmetric stresses has ℓ₇ = ℓ₈ = 0 identically — measured, not assumed.

𝔸_εε is dimensionless, so no scale divides out: see dimensionless_scale.

MeanFieldHomogenization.NeuralInclusions.StressLocTI Type

Same six-component transversely isotropic form as StrainLocTI, for the stress-side localization tensor 𝔸_σε.

The two differ only in their physical dimension, and therefore in what makes them dimensionless: 𝔸_εε is of degree 0 in the moduli, 𝔸_σε of degree +1.

MeanFieldHomogenization.NeuralInclusions.StrainLocCubic Type

Cubic strain localization tensor: the three components (α, β, γ) of a TensCubic on the projectors (𝕁, 𝔼, 𝕋), for a morphology left invariant by the octahedral group — a supersphere, a cubic array.

Three and not six, where StrainLocTI needs six for the same object one class up. A localization tensor has no major symmetry in general, which is exactly why the TI case cannot use the five-component major-symmetric Walpole form; but a tensor with the minor symmetries and cubic symmetry is major-symmetric automatically, the three projectors being symmetric. So nothing is dropped here: the class simply has no antisymmetric content to lose.

Dimensionless, like its TI sibling: 𝔸_εε is of degree 0 in the reference moduli.

MeanFieldHomogenization.NeuralInclusions.GradLocISO2 Type

Gradient localization tensor of a cavity in transport: the single component of an isotropic 2nd-order tensor, for gradient_gradient_loc.

One component, whatever the shape's symmetry class, because a second-order tensor invariant under the octahedral group is isotropic — the same fact that leaves the cubic class with no eigenstructure to key an orientation on. So a supersphere, which needs three constants in elasticity, needs one here, and the frame is irrelevant.

Distinct from HillISO2, which has the same one component and a different dimension. 𝑨_∇∇ is of degree 0 in the reference conductivity, so nothing divides out; the Hill tensor is of degree −1, so k₀ does. Reusing HillISO2 for a localization tensor would therefore divide the prediction by k₀ and be silently wrong at every k₀ ≠ 1 — see dimensionless_scale.

MeanFieldHomogenization.NeuralInclusions.GradLocTI2 Type

Gradient localization tensor of a cavity in transport, transversely isotropic: the two components (a, b) of TensTI{2} about the axis, for gradient_gradient_loc on an axisymmetric morphology.

Two where GradLocISO2 has one, because a body of revolution distinguishes its axis: a superspheroidal pore has R₁₁ ≠ R₃₃ for every aspect ratio but the sphere.

Distinct from HillTI2 for the same reason GradLocISO2 is distinct from HillISO2, and it is worth repeating because the two are otherwise interchangeable in shape: 𝑨_∇∇ is of degree 0 in the reference conductivity while the Hill tensor is of degree −1, so reusing HillTI2 would divide every prediction by k₀ and be silently wrong away from k₀ = 1.

The axis is column 3 of the inclusion basis, the package's usual convention and the one the axisymmetric solver revolves about.

The surrogate

MeanFieldHomogenization.NeuralInclusions.NeuralSurrogate Type

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.

MeanFieldHomogenization.NeuralInclusions.Provenance Type

Provenance

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

Field	Meaning
teacher	what generated the labels, e.g. "hill_tensor(Ellipsoid, TensISO)"
nsamples / nvalidation	training and held-out set sizes
max_block_error	worst error over the held-out set, relative to the tensor's own magnitude — the headline number
max_rel_error / rms_rel_error	per-component relative error, a diagnostic
epochs	how many passes the optimizer actually ran (after early stopping)
created	timestamp, for telling two retrainings apart
notes	free text

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.

MeanFieldHomogenization.NeuralInclusions.worst_error Function

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.

MeanFieldHomogenization.NeuralInclusions.check_domain Function

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.

MeanFieldHomogenization.NeuralInclusions.predict_components Function

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.

What the network predicts

The symmetry class, the major symmetry, the homogeneity in the reference moduli and the frame are enforced by these types rather than fitted — see What is exact, and what is fitted.

MeanFieldHomogenization.NeuralInclusions.AbstractHillClass Type

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.

MeanFieldHomogenization.NeuralInclusions.HillISO Type

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

MeanFieldHomogenization.NeuralInclusions.HillTI Type

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.

MeanFieldHomogenization.NeuralInclusions.HillOrtho Type

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.

MeanFieldHomogenization.NeuralInclusions.HillISO2 Type

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

MeanFieldHomogenization.NeuralInclusions.HillTI2 Type

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

MeanFieldHomogenization.NeuralInclusions.AbstractOutputSpec Type

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.

MeanFieldHomogenization.NeuralInclusions.DimensionlessHill Type

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.

MeanFieldHomogenization.NeuralInclusions.AnchoredHill Type

AnchoredHill(class, baseline::Symbol)

The network predicts the components of 𝕄 = 𝔸_b⁻¹ : 𝔸, where 𝔸_b is a closed-form baseline named by baseline. The decoder returns 𝔸_b : 𝕄, so the answer is exact wherever the baseline is exact, and the network only ever learns the departure from it.

This is the same idea as AffineHill — do not spend capacity on a dependence that is already known — applied to a different kind of knowledge. AffineHill factors out the material dependence exactly; AnchoredHill factors out a whole shape family on which a closed form exists.

When it pays, and when it does not — measured, so nobody repeats the experiment. The gain depends on how close the baseline is, not on its existence.

On the axisymmetric superspheroidal cavity it does not pay. A superspheroid at p = 1 is a spheroid, so an entire face of that box is exact — but the box reaches p = 0.25, which is far from it, and judged on the decoded tensor the anchored fit is worse: rms 1.30e-2 against 5.31e-3, median 9.16e-3 against 2.64e-3, gaining only 14 % on the extreme maximum. The reason is structural and worth knowing: near the exact face 𝕄 ≈ 𝕀, whose off-diagonal Walpole components are zero, so the anchored target has zeros by construction exactly where the anchor is perfect, and uniform relative effort in 𝕄 is not uniform in 𝔸.

Where it should pay is a baseline that is near over the whole box — a layered spheroid, whose homogeneous limit is reached as soon as the layers' moduli agree, so moderate contrast is a genuinely small departure and 𝕄 stays close to 𝕀 throughout.

The algebra stays in the Walpole basis, which is the whole reason this is a specification and not a division. 𝔸_b⁻¹ : 𝔸 is formed with inv and the double contraction on transversely isotropic tensors — a closed-form 2×2 inverse plus two scalars — where a component-wise ratio would mix the transverse block and manufacture sign changes that are an artifact of the reading.

Available baselines are the keys of anchor_baselines.

MeanFieldHomogenization.NeuralInclusions.anchor_baselines Function

anchor_baselines() -> Tuple{Symbol}

The baselines AnchoredHill can name. A baseline must be a closed form of the same symmetry class as the surrogate's, computable from the features alone plus the reference medium — otherwise it could not be evaluated at prediction time, when no solve is available.

• :spheroid_cavity — the exact cavity of revolution, 𝔸_b = (𝕀 − ℙ : ℂ₀)⁻¹, with the aspect ratio read off the features. Exact for a superspheroid at p = 1, and for a layered spheroid whose layers share the matrix's moduli.

MeanFieldHomogenization.NeuralInclusions.anchor_tensor Function

anchor_tensor(spec, x, features, P₀, frame) -> AbstractTens

The baseline tensor of an AnchoredHill specification at feature vector x. Evaluated from the features alone, so it is available at prediction time.

MeanFieldHomogenization.NeuralInclusions.encode Function

encode(spec, tensor, P₀, frame, x, features; atol) -> Vector

The inverse of decode: the components a label must carry for this specification. Used when building a dataset, and it is what keeps training and prediction in the same space — an anchored surrogate whose labels were the raw components would be silently wrong everywhere.

MeanFieldHomogenization.NeuralInclusions.spec_baseline Function

spec_baseline(spec) -> Union{Symbol, Nothing}

The anchor baseline's name, or nothing for a specification that has none. Serialized beside the class so a saved model round-trips.

MeanFieldHomogenization.NeuralInclusions.AffineHill Type

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.

MeanFieldHomogenization.NeuralInclusions.ncomponents Function

ncomponents(class) -> Int

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

MeanFieldHomogenization.NeuralInclusions.tensor_order Function

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.

MeanFieldHomogenization.NeuralInclusions.nterms Function

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.

MeanFieldHomogenization.NeuralInclusions.noutputs Function

noutputs(spec) -> Int

Width of the network's output layer.

MeanFieldHomogenization.NeuralInclusions.needs_nu Function

needs_nu(spec) -> Bool

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

MeanFieldHomogenization.NeuralInclusions.build Function

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.

MeanFieldHomogenization.NeuralInclusions.components Function

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.

MeanFieldHomogenization.NeuralInclusions.decode Function

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.

MeanFieldHomogenization.NeuralInclusions.material_coeffs Function

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.

MeanFieldHomogenization.NeuralInclusions.dimensionless_scale Function

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.

MeanFieldHomogenization.NeuralInclusions.hill_class Function

hill_class(name::Symbol) -> AbstractHillClass

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

MeanFieldHomogenization.NeuralInclusions.output_spec Function

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

Rebuild an output specification from its serialized names.

MeanFieldHomogenization.NeuralInclusions.apply_transform Function

apply_transform(kind, z)

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

MeanFieldHomogenization.NeuralInclusions.invert_transform Function

invert_transform(kind, y)

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

MeanFieldHomogenization.NeuralInclusions._feature Function

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

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

Name	Meaning	Range
:log_aspect	log(distinct axis / equal axes) of a spheroid	> 0 prolate, < 0 oblate, 0 sphere
:log_r2	log(a₂/a₁) of a sorted ellipsoid	≤ 0
:log_r32	log(a₃/a₂) of a sorted ellipsoid	≤ 0
:log_p	log of a supershape's concavity exponent	< 0 concave, 0 the ellipsoidal control
:nu0	Poisson ratio of the isotropic reference medium

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.

A SampleBox is linear, so the feature is the sampling law. That is the whole reason :log_p exists beside :p, which a pore also exposes. The resistivity contribution of an axisymmetric concave cavity runs from 1.66 at p = 0.6 to 15.0 at p = 0.20, near-diverging as the body tends to a crack pierced by a needle. Uniform in p, most of the sample budget lands where the response is flat and almost none where it turns over; uniform in log p it is spread evenly over the variation. The choice belongs to the study, not to the feature: :p remains available for a box whose response is well behaved.

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

MeanFieldHomogenization.NeuralInclusions.raw_features Function

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.

MeanFieldHomogenization.NeuralInclusions._class_frame Function

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

MeanFieldHomogenization.NeuralInclusions._canonical_axes Function

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

MeanFieldHomogenization.NeuralInclusions._spheroid_axis_index Function

_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°.

Sampling and labeling

MeanFieldHomogenization.NeuralInclusions.SampleBox Type

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

MeanFieldHomogenization.NeuralInclusions.Dataset Type

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.

MeanFieldHomogenization.NeuralInclusions.generate_dataset Function

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.

atol is the residual a label may leave when projected onto the class, and the default of 1.0e-8 is right for an analytic teacher, which lands in its class to round-off. A teacher that solves something does not: a finite-element cell measures a tensor carrying the discretization error of its mesh, and asking it for 1.0e-8 rejects a perfectly good label. Raise it to the accuracy the teacher actually has — for the octant cell at level 3 that is a few times 1.0e-5 — and keep it tight enough that a wrong class or a wrong frame, which misses by orders of magnitude more, still fails loudly.

reference overrides how the reference medium is built, as reference(box, x_full) -> AbstractTens. The default dispatches on the class, and for most classes that is the only sensible thing. It is not for StrainLocTI and StressLocTI: those describe the localization of a morphology, and whether the reference can be recovered from the features at all depends on the morphology and not on the class. A heterogeneous one carries its constituents inside itself, so scaling the reference changes the contrast and changes the answer — guessing there would train on corrupted labels, which is the failure this file exists to prevent. A cavity has no constituent and is of degree 0 in the reference, so the same class is perfectly usable, and the caller is the one who knows which case they are in. Hence a keyword rather than a method: it puts the statement where the knowledge is.

MeanFieldHomogenization.NeuralInclusions.sample_box Function

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.

MeanFieldHomogenization.NeuralInclusions.grid_box Function

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.

MeanFieldHomogenization.NeuralInclusions.halton Function

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

MeanFieldHomogenization.NeuralInclusions.feature_index Function

feature_index(box, name) -> Int

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

MeanFieldHomogenization.NeuralInclusions.fit_scaling Function

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.

MeanFieldHomogenization.NeuralInclusions.validate_surrogate Function

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.

MeanFieldHomogenization.NeuralInclusions.report_surrogate Function

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.

MeanFieldHomogenization.NeuralInclusions.component_labels Function

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 𝕍; for AnchoredHill they are the class's, the anchored target living in the same basis.

Training

train_surrogate is the seam of the MeanFieldHomogenizationLuxExt extension: the method below is the fallback that raises when the extension is not loaded.

MeanFieldHomogenization.NeuralInclusions.TrainingOptions Type

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

Option	Default	Meaning
hidden	[32, 32]	widths of the hidden layers; input and output widths follow from the feature list and the output specification
activation	:tanh	hidden activation — must be smooth, see mlp.jl
epochs	4000	maximum passes over the training set
batchsize	128	mini-batch size; 0 means full batch
learning_rate	1.0e-2	initial Adam step
decay	0.3	multiplicative decay applied to the step at each plateau
patience	250	epochs without validation improvement before decaying, and 3·patience before stopping
seed	20260730	RNG seed of the weight initialization, so a retraining is reproducible
verbose	true	print the validation curve as it goes

MeanFieldHomogenization.NeuralInclusions.train_surrogate Function

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.

MeanFieldHomogenization.NeuralInclusions.assemble_surrogate Function

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.

MeanFieldHomogenization.NeuralInclusions.network_widths Function

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.

The network, and its serialization

MeanFieldHomogenization.NeuralInclusions.MLP Type

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.

MeanFieldHomogenization.NeuralInclusions.NNDense Type

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.

MeanFieldHomogenization.NeuralInclusions.glorot_mlp Function

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

Fresh network with Glorot-uniform weights and zero biases: layer k draws from , 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.

MeanFieldHomogenization.NeuralInclusions.softplus Function

softplus(x)

Smooth positive activation  , evaluated in the numerically stable form    so that large |x| neither overflows nor loses the linear branch.

MeanFieldHomogenization.NeuralInclusions.activation Function

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.

MeanFieldHomogenization.NeuralInclusions.activation_name Function

activation_name(σ) -> Symbol

Reverse lookup, for serialization.

MeanFieldHomogenization.NeuralInclusions.layer_widths Function

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.

MeanFieldHomogenization.NeuralInclusions.layer_activations Function

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

Activation name of each layer, output layer included.

MeanFieldHomogenization.NeuralInclusions.nparams Function

nparams(m::MLP) -> Int

Total number of weights and biases.

MeanFieldHomogenization.NeuralInclusions.save_surrogate Function

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.

MeanFieldHomogenization.NeuralInclusions.load_surrogate Function

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.

MeanFieldHomogenization.NeuralInclusions.model_path Function

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.

MeanFieldHomogenization.NeuralInclusions.shipped_models Function

shipped_models() -> Vector{String}

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

MeanFieldHomogenization.NeuralInclusions.SURROGATE_FORMAT Constant

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.

MeanFieldHomogenization.NeuralInclusions.MODEL_DIR Constant

MODEL_DIR

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