API — Schemes

MeanFieldHomogenization.Schemes Module

MeanFieldHomogenization.Schemes

Mean-field homogenization schemes. Provides the RVE container (matrix + named phases with their volume fractions or crack densities, plus an optional distribution shape) and the suite of homogenization HomogenizationScheme types: bounds (Voigt, Reuss), one-shot schemes with a matrix (Dilute, DiluteDual, MoriTanaka, Maxwell, PonteCastanedaWillis), iterative self-consistent schemes (SelfConsistent, AsymmetricSelfConsistent) and the differential scheme (Differential) with user-selectable trajectory.

Public entry point: homogenize(rve, scheme; property=:C). The scheme can also be passed as a Symbol shortcut (:MT, :SC, …).

RVE / Phase / Amount

MeanFieldHomogenization.Schemes.RVE Type

RVE{T<:Number, S<:Union{Nothing,AbstractDistributionShape}}

Multi-phase representative volume element — a description of a microstructure, and nothing more. No phase is singled out: which one, if any, plays the reference medium of a homogenization scheme is stated on that scheme (see matrix_name), because the same RVE is a matrix/inclusion composite under Mori-Tanaka and a matrix-free aggregate under the self-consistent scheme.

Fields:

• phase_names::Vector{Symbol} — phases in declaration order.

• phases::Dict{Symbol,Phase} — geometry + properties of each phase.

• amounts::Dict{Symbol,AbstractAmount} — VolumeFraction, CrackDensity or Remainder of each phase. Each entry keeps its own element type.

• distribution_shape::S — outer envelope used by Maxwell / PCW, wrapped in UniformDistribution; nothing when the caller declared none, which is the default. Only those two schemes read it, the way only some schemes read a phase's geometry — but unlike a geometry it has no harmless default, so they raise rather than assume one (see distribution_shape).

• closure — the AbstractFractionClosure turning the declared fractions into the ones the schemes use; nothing until inferred.

• rest_name — the phase declared fraction = :rest, if any.

• f_sum — cached Σ f, maintained by add_phase! and by set_amount!. Do not write amounts directly: that would leave the cache stale.

T is the declared amount element type, i.e. a floor for promotion, not a constraint: it seeds zero/one for an RVE whose amounts are all narrower (an Int fraction still yields a Float64 matrix fraction under the default T = Float64), and any amount handed to add_phase! that is wider than T is stored as such rather than converted down. Amounts, moduli and geometries therefore each carry their own element type, promoted only where the values meet:

rve = RVE()                                          # nothing to declare
add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => C_complex); fraction = :rest)    # complex moduli
add_phase!(rve, :I, Ellipsoid(1.0), Dict(:C => C_complex);
           fraction = 0.3)                             # real fraction
add_phase!(rve, :J, Ellipsoid(1.0), Dict(:C => C1);
           fraction = dual_x)                          # Dual fraction

eltype(rve) reports the effective element type (the promotion of the floor with every stored amount); eltype(RVE{T,S}) reports the declared floor T.

Construction is two-step:

rve = RVE(; distribution_shape = nothing)   # or RVE{ComplexF64}()
add_phase!(rve, :M, ellipsoid_matrix, Dict(:C => C0); fraction = :rest)
add_phase!(rve, :I1, ellipsoid_inc, Dict(:C => C1); fraction = 0.2)
add_phase!(rve, :CRACK, penny_crack, Dict(:C => C0); density = 0.05)

See also add_phase!, volume_fraction, AbstractFractionClosure, validate_rve.

MeanFieldHomogenization.Schemes.Phase Type

Phase(geometry::AbstractInclusion, properties::Dict{Symbol,<:AbstractTens})

A single phase of a RVE: one inclusion geometry (ellipsoid, crack, …) together with one or several material property tensors indexed by symbol (:C for stiffness, :K for conductivity, …).

The geometry is field-typed AbstractInclusion (rather than parametric Phase{I}) so that a heterogeneous RVE mixing ellipsoids and cracks can be stored in a single Dict{Symbol,Phase} without losing information at construction time. Specialization happens at the dispatch site (hill_tensor(phase.geometry, …), cod_tensor(phase.geometry, …)).

MeanFieldHomogenization.Schemes.AbstractAmount Type

AbstractAmount{T<:Number}

Supertype for the quantity attached to a phase in a RVE. Two concrete subtypes:

• VolumeFraction — for solid (ellipsoidal) inclusions and the matrix; obeys the unit-sum constraint f_matrix = 1 - Σ_other VolumeFraction.value.

• CrackDensity — for flat cracks (Budiansky-O'Connell density); does not participate in the unit-sum constraint, since the volume contribution of a flat crack vanishes in the penny limit while the density remains finite.

The type parameter T is the element type of the stored value (Float64, ForwardDiff.Dual, Complex{Float64}, …) and propagates through every scheme that consumes the amount. Each amount carries its own T: a RVE may hold a Float64 fraction next to a ForwardDiff.Dual one, the element types being promoted where the values actually meet (see RVE).

MeanFieldHomogenization.Schemes.VolumeFraction Type

VolumeFraction(f) <: AbstractAmount

Volume fraction of a solid inclusion (or of the matrix).

MeanFieldHomogenization.Schemes.CrackDensity Type

CrackDensity(ε) <: AbstractAmount

Budiansky-O'Connell crack density of a population of flat cracks.

MeanFieldHomogenization.Schemes.AbstractDistributionShape Type

AbstractDistributionShape

Supertype for the outer envelope of the phase distribution used by the Maxwell and PonteCastanedaWillis schemes.

It is a field of the RVE rather than of those schemes because it is microstructure: the ellipsoidal symmetry of the medium's two-point statistics ([28]), a measurable property like a phase's shape. That is what distinguishes it from the reference medium, which is a modeling decision and therefore lives on the scheme. It has no default, though — see distribution_shape.

Currently a single concrete subtype is shipped:

• UniformDistribution — a single shape applied to every inclusion phase (Maxwell 1873 ; Ponte-Castañeda & Willis 1995).

Future extension (placeholder, not implemented in this PR): a PairwiseDistribution carrying a per-pair (i, j) ↦ shape mapping ([29]). Adding it will only require a new concrete subtype + matching _evaluate(rve, ::Maxwell|::PonteCastanedaWillis, …) methods — no public-API change.

MeanFieldHomogenization.Schemes.UniformDistribution Type

UniformDistribution(shape::AbstractInclusion) <: AbstractDistributionShape

Single distribution shape applied to every inclusion phase. The argument-free constructor UniformDistribution() returns a unit sphere — an isotropic distribution, which makes Maxwell and PonteCastanedaWillis coincide exactly with Mori-Tanaka. That is a legitimate choice and has to be made explicitly: an RVE given no distribution shape keeps nothing and those two schemes raise rather than assume a sphere.

MeanFieldHomogenization.Schemes.AbstractSymmetrize Type

AbstractSymmetrize

Specifies how a phase's localization tensor (and the derived stiffness / compliance / conductivity / resistivity contributions) is averaged over an orientation distribution before being used in the homogenization formula.

Three concrete subtypes are shipped :

• NoSymmetrize (default) — keep the contribution as computed for the single-orientation inclusion stored in the phase.

• IsoSymmetrize — average over all rotations (uniform spatial distribution of orientations) ; produces an isotropic projection.

• TISymmetrize — average over rotations around a specified axis (uniaxial uniform distribution) ; produces a transversely-isotropic projection.

This mirrors C++ ECHOES's symmetrize=[ISO] / symmetrize=[TI] keyword on ellipsoid(), but moved to the RVE side (just like volume fractions) : the same inclusion type can be re-used in different RVEs with different distribution assumptions, and a single inclusion remains usable for localization-tensor calculations without any RVE.

MeanFieldHomogenization.Schemes.NoSymmetrize Type

NoSymmetrize() <: AbstractSymmetrize

Default. The localization tensor is used as computed for the single orientation defined by the inclusion's basis.

MeanFieldHomogenization.Schemes.IsoSymmetrize Type

IsoSymmetrize() <: AbstractSymmetrize

The localization tensor is averaged over a uniform spatial distribution of orientations, equivalent to projecting onto the isotropic basis (J, K_proj) for 4th-order tensors and onto the spherical part for 2nd-order tensors. Produces an isotropic phase contribution regardless of the inclusion's actual shape.

MeanFieldHomogenization.Schemes.TISymmetrize Type

TISymmetrize(axis = (0, 0, 1); reference_projection = :iso) <: AbstractSymmetrize

The localization tensor is averaged exactly over rotations about axis (uniaxial uniform distribution). The average preserves the full axially-invariant structure — including the non-major-symmetric components of concentration tensors — and returns a TensND.TensTI{4,T,8} (resp. TensTI{2,T,3} at 2nd order).

reference_projection controls how the reference medium is pre-projected before the localization tensor of this phase is computed :

• :iso (default) — project the reference to its isotropic average. Every inclusion orientation then has an analytical, ForwardDiff-compatible Hill branch. This is an approximation whenever the reference is not isotropic; it is exact at the isotropic fixed point of a self-consistent iteration.

• :none — use the reference as is. Non-coaxial anisotropic references route through the general-anisotropy Hill branch (NestedQuadGK, ForwardDiff-compatible but quadrature-priced).

• :ti — project the reference to its best-fit TI form about axis. Only valid when the phase's inclusions are coaxial with axis (the TI-coaxial analytical Hill branch applies).

MeanFieldHomogenization.Schemes.phase_symmetrize Function

phase_symmetrize(rve, name::Symbol) -> AbstractSymmetrize

Return the orientation-distribution projection declared for phase name. Defaults to NoSymmetrize if none was set.

MeanFieldHomogenization.Schemes.add_phase! Function

add_phase!(rve, name::Symbol, geometry, properties::AbstractDict;
           fraction = nothing, density = nothing, symmetrize = nothing)

Register an inclusion phase with the given geometry and material properties. Exactly one of fraction (for ellipsoidal inclusions and solid inhomogeneities) or density (for cracks) must be supplied; fraction produces a VolumeFraction, density a CrackDensity.

Both fraction and density are stored under the RVE's declared element-type floor T: widened to promote_type(T, typeof(value)), never narrowed. A complex, ForwardDiff.Dual or symbolic amount is therefore accepted by a plain RVE(), and phases may carry amounts of different element types.

The optional symmetrize kwarg declares an orientation-distribution projection of this phase's localization tensor : :iso (uniform spatial distribution → isotropic projection), :ti (uniaxial uniform around z-axis), TISymmetrize(axis) (around an arbitrary axis), or pass an explicit AbstractSymmetrize instance. The default NoSymmetrize keeps the inclusion's actual single-orientation tensor.

MeanFieldHomogenization.Schemes.phase_names Function

phase_names(rve::RVE) -> Vector{Symbol}

Names of every phase, in declaration order. This is what a scheme that distinguishes no phase — the self-consistent family, the bounds — iterates over; schemes built on a reference medium use inclusion_phase_names instead.

MeanFieldHomogenization.Schemes.inclusion_phase_names Function

inclusion_phase_names(rve::RVE, matrix::Symbol) -> Vector{Symbol}

Names of the phases other than matrix, in declaration order — the phases a matrix-based scheme localizes in its reference medium.

The phase to exclude has to be named: the reference medium is a property of the scheme, and an RVE holds no such designation of its own. Use phase_names for the schemes that distinguish none.

MeanFieldHomogenization.Schemes.phase_property Function

phase_property(rve, name::Symbol, key::Symbol) -> AbstractTens

Return the property tensor key (e.g. :C, :K) of the phase named name.

If the stored value is a Homogenized — a declaratively nested cell — it is resolved here, memoized for the duration of the enclosing homogenize call, so that every scheme kernel sees a plain tensor without knowing that nesting exists. Use phase_property_raw to inspect the stored value without triggering that resolution.

MeanFieldHomogenization.Schemes.volume_fraction Function

volume_fraction(rve, name::Symbol) -> Number

Volume fraction of phase name, as the schemes see it: the RVE's AbstractFractionClosure has already been applied, so a Remainder reads as 1 - Σ f and, under RescaledFractions, a declared fraction reads as its renormalized value. This is the one accessor a scheme should use.

The declared value — what a parameter lens round-trips — is get_param(rve, AmountParameter(name)) instead. The two differ under RescaledFractions, and only there.

A CrackDensity phase reads as zero: a flat crack carries no volume.

MeanFieldHomogenization.Schemes.crack_density Function

crack_density(rve, name::Symbol) -> Number

Crack density of phase name. Returns a zero of the phase's own element type (promoted with the RVE floor) if the phase carries a VolumeFraction instead of a CrackDensity.

MeanFieldHomogenization.Schemes.remainder_phase_name Function

remainder_phase_name(rve::RVE) -> Union{Nothing, Symbol}

The phase declared fraction = :rest, or nothing when the RVE designates none.

MeanFieldHomogenization.Schemes.remainder_volume_fraction Function

remainder_volume_fraction(rve::RVE) -> Number

The volume left over by the declared fractions, 1 - Σ f — what the phase declared fraction = :rest takes up.

CrackDensity entries do not contribute: a flat crack's volume vanishes in the penny limit while its density stays finite.

This is a read of the cached f_sum field plus one subtraction, not a recomputation. The unit comes from the accumulator rather than from the declared floor, so a symbolic RVE yields 1 - f and not 1.0 - f.

MeanFieldHomogenization.Schemes.phase_property_raw Function

phase_property_raw(rve, name::Symbol, key::Symbol)

The value stored under key on phase name, without resolving a Homogenized. Type inspections must use this — resolving would run a whole inner homogenization just to look at a type.

MeanFieldHomogenization.Schemes.validate_rve Function

validate_rve(rve::RVE)

Sanity-check the RVE: at least one phase, no negative amount, and the fraction-closure policy satisfied. Throws ArgumentError on hard failures; emits @warn for a negative complement under the default ComplementFraction(on_negative = :warn) — a non-physical RVE stays useful for symbolic or Dual exploration, but is flagged.

MeanFieldHomogenization.Schemes.promote_rve Function

promote_rve(rve, ::Type{T}) -> RVE

Return a copy of rve whose declared floor is T and whose amounts are all converted to promote_type(T, ·). Rarely needed — amounts promote themselves where they are consumed — but useful to force an element type on an RVE built elsewhere.

Fraction closures

MeanFieldHomogenization.Schemes.Remainder Type

Remainder() <: AbstractAmount

Marker amount of the phase whose volume fraction is derived rather than declared — the one that takes up whatever the others leave.

It carries no value: under ComplementFraction the value is 1 - Σ f, a property of the RVE and not of the phase. Read the resolved number with volume_fraction; amount_value deliberately throws, so a loop that assumed every amount carries a number fails by name instead of silently skipping the phase.

Its element type is Union{}, which promotes away: an RVE holding a Remainder next to Float64 and ForwardDiff.Dual fractions has exactly the eltype it would have without it.

MeanFieldHomogenization.Schemes.AbstractFractionClosure Type

AbstractFractionClosure

How an RVE turns the volume fractions the caller declared into the fractions the schemes use. Three concrete policies ship:

• ComplementFraction — one phase, declared fraction = :rest, absorbs 1 - Σ f;

• RescaledFractions — every declared fraction is divided by their sum, so relative proportions may be given;

• StrictFractions — the declared fractions must already sum to one.

The policy is inferred when it is not given: ComplementFraction() as soon as a phase is declared with fraction = :rest, StrictFractions() otherwise.

CrackDensity never takes part. A flat crack has no volume, so it is outside the unit-sum constraint, is never renormalized, and is never absorbed into a complement.

MeanFieldHomogenization.Schemes.StrictFractions Type

StrictFractions(; atol = 1e-10) <: AbstractFractionClosure

The declared volume fractions must sum to one within atol; anything else is an error naming the sum. The default policy for an RVE in which no phase is declared with fraction = :rest — a polycrystal, a granular aggregate, any microstructure in which no phase is a leftover.

MeanFieldHomogenization.Schemes.ComplementFraction Type

ComplementFraction(; on_negative = :warn) <: AbstractFractionClosure

The single phase declared fraction = :rest takes the volume the others leave, 1 - Σ f.

on_negative says what happens when that complement comes out negative, i.e. when the declared inclusion fractions already exceed one: :warn (default, the historical behavior — a non-physical RVE is still useful for symbolic or ForwardDiff exploration) or :error.

MeanFieldHomogenization.Schemes.RescaledFractions Type

RescaledFractions() <: AbstractFractionClosure

Every declared volume fraction is divided by their sum, so that only their ratios matter: 2, 3, 5 and 0.2, 0.3, 0.5 describe the same RVE.

No phase may then be declared fraction = :rest — there is no complement left to absorb. The rescaling is a plain division, so it differentiates: under this policy ∂C/∂f_i is the derivative along the normalized simplex, and raising one fraction lowers the others.

What a scheme requires of the RVE

The reference medium is named on the scheme; the distribution shape is declared on the RVE. Both are resolved before any kernel runs, and neither is guessed.

MeanFieldHomogenization.Schemes.matrix_name Function

matrix_name(scheme, rve) -> Symbol

The phase of rve that plays the reference medium of scheme.

Named explicitly on the scheme when the caller said so (MoriTanaka(:cem)); otherwise the phase that absorbs the volume complement, when the RVE designates one. Choosing a matrix is a modeling decision — a different choice is a different composite, not a rounding difference — so an RVE that designates none and a scheme that names none is an error, and the message lists the candidates rather than picking for the caller.

MeanFieldHomogenization.Schemes.reference_property Function

reference_property(rve, scheme, key::Symbol) -> AbstractTens

Property key of the phase that plays scheme's reference medium — the scheme-aware replacement for reading a property off "the matrix" of an RVE.

MeanFieldHomogenization.Schemes.requires_matrix Function

requires_matrix(scheme) -> Bool

Whether scheme computes localization tensors in a reference medium taken from one distinguished phase, and therefore cannot be evaluated until that phase is decided.

false for the bounds, for the self-consistent scheme and for Laminated; true for the matrix-based estimates.

MeanFieldHomogenization.Schemes.scheme_matrix Function

scheme_matrix(scheme) -> Union{Nothing, Symbol}

The phase the caller named as the reference medium of scheme, or nothing when the scheme has no such notion or the caller left the choice open.

Generic fallback: nothing. A scheme carrying a matrix field overrides it.

MeanFieldHomogenization.Schemes.host_phase_name Function

host_phase_name(rve, declared, who::AbstractString) -> Symbol

The phase that plays the continuous host of rve for a consumer that is not a homogenization scheme — the solid skeleton of a fracture network, the homogeneous solid of micro-poromechanics.

declared is what the caller named, or nothing to fall back on the phase that absorbs the volume complement. who names the caller in the error message: an RVE that designates no complement leaves the host undecided, and which phase is the skeleton is a modeling statement no container can make on the caller's behalf.

MeanFieldHomogenization.Schemes.distribution_shape Function

distribution_shape(rve, scheme) -> AbstractDistributionShape

The outer envelope of the spatial distribution that scheme reads off rve.

An RVE that declares none is an error, and deliberately so. A spherical distribution — the value this used to default to — makes both Maxwell and PonteCastanedaWillis reduce exactly to Mori-Tanaka. Supplying one silently therefore answered a scheme whose whole purpose is a non-spherical distribution with the estimate it generalizes, and did so without a word. A sphere remains a perfectly legitimate choice; it just has to be the caller's, not the container's.

MeanFieldHomogenization.Schemes.requires_distribution_shape Function

requires_distribution_shape(scheme) -> Bool

Whether scheme builds its estimate on the shape of the spatial distribution of the inclusions, and therefore cannot be evaluated until the RVE declares one.

true for Maxwell and PonteCastanedaWillis; false for every other scheme, which ignores the field the way Voigt ignores a phase's geometry.

The self-consistent solvers

Internal, but the normative statement of what abstol and reltol ask for — see Solver tolerances for the prose.

MeanFieldHomogenization.Schemes._solve_sc Function

_solve_sc(algo, step, x0; abstol, reltol, maxiters, damping, verbose,
          select_best, kw...) -> AbstractTens

Generic solver dispatcher for SC fixed points. Built-in:

• AndersonDefault — Picard with relaxation (x_{n+1} = (1-damping)·step(x_n) + damping·x_n). damping = 0.0 default; raise to ≈ 0.5 for high-contrast iterations that overshoot. Convergence near a bifurcation (e.g. SC at the porous-percolation threshold) is intrinsically slow because the Picard Jacobian eigenvalue approaches 1 there; in that regime, set select_best = true to return the best iterate observed during the loop, or load NonlinearSolve.jl and switch to Newton/Anderson via the algorithm keyword.

• NewtonDefault — built-in Newton-Raphson with a ForwardDiff Jacobian on the canonical symmetry components and an Armijo line search. Dependency-free; ships with the package.

• AutoNonlinear — auto-resolving marker: uses a globalized NonlinearSolve.jl algorithm (TrustRegion) when the weak extension MeanFieldHomogenizationNonlinearSolveExt is loaded (using NonlinearSolve), else falls back to NewtonDefault.

Any algorithm from NonlinearSolve.jl (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), …) can also be passed directly as algorithm; it is handled by the same weak extension, through a ForwardDiff-safe implicit-function-theorem lift so that derivative/gradient/jacobian (see sensitivities.jl) work transparently regardless of which solver is selected.

Convergence is declared when ‖x_new − x_old‖ ≤ abstol + reltol · ‖x_old‖ (absolute and relative tolerance, additive convention; pass abstol = 0 to require purely relative convergence). Default values: abstol = 1e-12, reltol = 1e-8.

‖·‖ is the Frobenius norm of the tensor, for every solver, so that a given abstol expresses one requirement regardless of algorithm and of the symmetry class the fixed point lives in. The Picard loop measures it directly; the Newton path reaches it by working in the isometric parametrization of _sc_param_weights; the NonlinearSolve path passes those same weights to SciML as the internalnorm of its termination condition, leaving the unknowns alone — a trust region is not invariant under a rescaling of them, and moving its metric moves where it stops.

Because a stiffness carries a physical magnitude, reltol · ‖x‖ dominates the sum at the defaults above: tightening abstol alone does not tighten the iteration. Set both, or abstol = 0, when a converged value is read off rather than plotted.

When select_best = true, the solver tracks the best iterate seen during the loop (smallest residual on the value field) and returns it at the end. Useful for high-contrast iterations where Picard oscillates around the fixed point: the last iterate may be worse than an earlier one. Default is false (return last iterate).

Non-convergence is reported via @debug (silent by default; set JULIA_DEBUG=MeanFieldHomogenization to surface it) rather than @warn. Near bifurcation points the Picard step intrinsically slows down (the linearized step has a Jacobian eigenvalue ≈ 1) and the residual stalls above tol_eff while still being negligibly small compared to the matrix-property scale; the returned iterate is informative even when the strict tolerance is not reached, so a default warning would be noise.

_solve_sc(::NewtonDefault, step, x0::AbstractTens; …) -> AbstractTens

Built-in Newton-Raphson SC solver, parameterizing the iterating estimate by its symmetry-class canonical components (TensND.get_data → (α, β) for iso, (ℓ₁, …, ℓ₆) for TI / Walpole, 9 components for ortho). At each Newton step:

1. Build the residual F(p) = canonical(step(rebuild(p))) − p,

2. Compute the Jacobian J = ∂F/∂p via ForwardDiff.jacobian,

3. Take the Newton step Δp = −J⁻¹·F(p) with backtracking line search (Armijo with shrinking factor 1/2, minimum step 1e-6).

4. Fall back to a single Picard step when the line search fails.

Compared to the SciML weak-extension path, this is dependency-free and specialized to the small parameter spaces of MeanFieldHomogenization symmetry classes (≤ 21 components for the most general aniso 4-tensor); the Jacobian is computed once per iteration through the same step function the AndersonDefault loop calls.

_solve_sc(::AutoNonlinear, step, x0::AbstractTens; kw...) -> AbstractTens

Resolver for AutoNonlinear: checks at runtime whether the weak extension MeanFieldHomogenizationNonlinearSolveExt is loaded (Base.get_extension) and, if so, delegates to ext.default_solve_sc(step, x0; kw...) — a globalized SciML algorithm (NonlinearSolve.TrustRegion()) run through the same ForwardDiff-safe (implicit-function-theorem) path as any other NonlinearSolve.jl algorithm. Falls back to the built-in NewtonDefault solver when the extension is not loaded, so SelfConsistent(algorithm = AutoNonlinear()) always works — with or without using NonlinearSolve.

MeanFieldHomogenization.Schemes._sc_param_weights Function

_sc_param_weights(prototype) -> (w::Vector{Float64}, isometric::Bool)

Per-component weights making the Euclidean norm of the Newton / SciML parameter vector equal the Frobenius norm of the tensor it encodes, so that abstol and reltol express the same requirement whatever the algorithm and whatever the symmetry class.

The canonical components are coordinates in a basis of the symmetry class ((α, β) on (𝕁, 𝕂) for TensISO, the Walpole coefficients for TensTI, …). Those bases are orthogonal but not orthonormal — ‖𝕁‖_F = 1 while ‖𝕂‖_F = √5 — so the plain Euclidean norm of the component vector understates the tensor norm, and understates it by a class-dependent factor. Scaling component i by wᵢ = ‖rebuild(eᵢ)‖_F restores the isometry exactly:

Newton is invariant under a diagonal rescaling of its unknowns — J becomes W J W⁻¹, the step W δ, and the iterate W(p + δ) — so the root and the sequence of iterates are untouched; only the norm the stopping test and the Armijo condition are measured in changes. A trust-region method is not invariant, and there the rescaling additionally puts its region in a physically meaningful metric.

The isometry needs the class basis to be orthogonal, which is checked rather than assumed: isometric = false (with unit weights) is returned when the probe fails — the general Mandel66 fallback among others — and the caller then measures the tensor norm directly instead.

Schemes

MeanFieldHomogenization.Schemes.HomogenizationScheme Type

HomogenizationScheme

Supertype for every mean-field homogenization scheme. Concrete subtypes:

• bounds — Voigt, Reuss ;

• one-shot with matrix — Dilute, DiluteDual, MoriTanaka, Maxwell, PonteCastanedaWillis ;

• iterative — SelfConsistent, AsymmetricSelfConsistent ;

• trajectory-based — DifferentialScheme.

MeanFieldHomogenization.Schemes.Voigt Type

Voigt() <: HomogenizationScheme

Voigt (uniform-strain) upper bound: .

MeanFieldHomogenization.Schemes.Reuss Type

Reuss() <: HomogenizationScheme

Reuss (uniform-stress) lower bound: .

MeanFieldHomogenization.Schemes.Laminated Type

Laminated() <: HomogenizationScheme

Exact solution of the periodic multilayer unit cell — a stack of parallel layers of common normal n, with no matrix and no reference medium. Applies to a Laminate cell, not to an RVE.

Unlike every other scheme in this file it is not an estimate: for a laminate it is the answer,

with ℙ_i the flat-inclusion Hill tensor of layer i and ℚ_i = ℂ_i − ℂ_i:ℙ_i:ℂ_i. Serves elasticity and transport, and accepts imperfect interfaces of spring / membrane / Kapitza / surface-conductive type. See the theory page on the laminate and [61].

MeanFieldHomogenization.Schemes.Dilute Type

Dilute() <: HomogenizationScheme

Dilute scheme:    where     is the size-independent stiffness contribution ([6]; [24]).

MeanFieldHomogenization.Schemes.DiluteDual Type

DiluteDual() <: HomogenizationScheme

Dual dilute scheme on the compliance:   , returning  .

MeanFieldHomogenization.Schemes.MoriTanaka Type

MoriTanaka() <: HomogenizationScheme

Mori-Tanaka scheme ([26]; [27]).

MeanFieldHomogenization.Schemes.Maxwell Type

Maxwell() <: HomogenizationScheme

Maxwell homogenization, using the RVE's distribution shape as the reference for the Hill polarization tensor.

MeanFieldHomogenization.Schemes.PonteCastanedaWillis Type

PonteCastanedaWillis() <: HomogenizationScheme

Ponte-Castañeda & Willis 1995 scheme — distribution-shape-aware generalization of Mori-Tanaka.

MeanFieldHomogenization.Schemes.SelfConsistent Type

SelfConsistent(; algorithm = AndersonDefault(), init = nothing, kwargs...) <: HomogenizationScheme

Self-consistent scheme — the morphology in which no phase is a matrix, so this scheme carries no reference-medium field at all and runs on an RVE that designates none.

The algorithm selects the non-linear solver; default is the built-in Anderson acceleration. Any solver from the SciML NonlinearSolve.jl package can be passed once using NonlinearSolve activates the weak extension MeanFieldHomogenizationNonlinearSolveExt.

init is the starting iterate. Unset, it is the phase taking up the volume complement when the RVE names one, and the Voigt average otherwise — the latter being what lets an RVE with no complement phase be solved at all. It also accepts :voigt, :reuss, a phase name, or an explicit tensor. A phase named :voigt or :reuss wins over the keyword: a phase name is the caller's data.

The seed picks the path, not the answer

Away from percolation the fixed point is unique and every seed reaches it: measured over a porous RVE at f ∈ {0.1, 0.3, 0.45, 0.6} and a stiff one at f ∈ {0.1, 0.3, 0.5, 0.7}, a Voigt seed and a phase seed agree to 1e-9 or better. They part only next to the porous percolation threshold (f = 0.49, 0.52), where the stiffness has collapsed to ~1e-4 against a solid at 72 and the damped Picard crawls along a degenerate fixed point; NewtonDefault agrees to 1e-7 even there. On a deeply percolated oblate RVE the Voigt start also makes TrustRegion stop where the finite-difference IFT Jacobian is exactly singular. Hence the default above.

Standard kwargs forwarded to the solver: abstol, reltol, maxiters, damping, verbose, select_best. Convergence is declared when ‖Δx‖ ≤ abstol + reltol · ‖x‖ (defaults abstol = 1e-12, reltol = 1e-8); pass abstol = 0 for a purely relative test. Because the stiffness carries a physical magnitude, it is reltol and not abstol that binds in the usual case — see _solve_sc and the schemes manual.

MeanFieldHomogenization.Schemes.AsymmetricSelfConsistent Type

AsymmetricSelfConsistent(; algorithm = AndersonDefault(), kwargs...) <: HomogenizationScheme

Asymmetric self-consistent scheme: iterates in stiffness or compliance space depending on the matrix-vs-Voigt-bound contrast, providing a better behavior than SelfConsistent in matrix-stiff / inclusion-soft regimes.

algorithm and init mean exactly what they do for SelfConsistent, and the same solver kwargs apply — abstol, reltol, maxiters, damping, verbose, select_best, with the same ‖Δx‖ ≤ abstol + reltol · ‖x‖ test and the same defaults. Unlike SelfConsistent it does carry a reference medium, because deciding between the stiffness and the compliance form is a comparison against that phase.

MeanFieldHomogenization.Schemes.DifferentialScheme Type

DifferentialScheme(; trajectory = Proportional(), nsteps::Int = 100,
                     abstol::Real = 1e-8, reltol::Real = 1e-6,
                     alg = nothing, formulation = :stiffness, kwargs...)

Differential scheme : integrates the Norris ODE on the fictitious incorporation time τ ∈ [0, 1] ([31]) :

with the volume balance df = (𝟙 - f ⊗ 𝐔)·dφ inverted by Sherman- Morrison so the user supplies effective volume fractions f_α(τ) along the chosen trajectory.

Keyword arguments

• trajectory — one of Proportional, Sequential, CustomPath, Path. Default Proportional().

• formulation — :stiffness (default) integrates the ODE above; :compliance integrates its exact dual    and inverts the result, so both return the same declared property. The two agree analytically (ℍ = −𝕊 : 𝐍 : 𝕊) and differ only in which variable carries the solver's error control: prefer :compliance for a medium softening towards percolation (porous, cracked), :stiffness for a stiffening one.

• nsteps — density of save points along τ (passed as saveat to the SciML ODE solver). The integration step is controlled by abstol / reltol, not by nsteps. See differential_path to read the saved states back.

• abstol, reltol — ODE solver tolerances (forwarded to OrdinaryDiffEq.solve).

• alg — explicit ODE algorithm. nothing selects Tsit5() (5th order adaptive Runge-Kutta). Pass any OrdinaryDiffEqAlgorithm instance to override (e.g. Vern9() for higher accuracy). Implicit algorithms must be built with a non-AD Jacobian (Rosenbrock23(autodiff = AutoFiniteDiff())): the RHS calls the Hill-tensor backends, which are not differentiable with respect to the ODE state.

• kwargs... — any other keyword is forwarded verbatim to OrdinaryDiffEq.solve (maxiters, dtmax, dt, callback, …).

MeanFieldHomogenization.Schemes.DifferentialTrajectory Type

DifferentialTrajectory

Supertype describing the path through the multi-phase volume-fraction space used by the DifferentialScheme scheme. Concrete subtypes:

• Proportional (default) — every phase grows linearly with the fictitious incorporation time τ ∈ [0, 1], all phases reaching their target simultaneously.

• Sequential — phases are introduced one after the other in the user-supplied order, each occupying a contiguous slice of τ.

• CustomPath — explicit per-phase trajectory as a vector of monotone non-decreasing values; piecewise-linear interpolated along τ ∈ [0, 1].

• Path — explicit per-phase trajectory as a callable τ -> f(τ) (auto-differentiated by ForwardDiff); the natural API for the multi-phase incorporation-sequence ODE ([31]; the user's hand-written DEM note).

MeanFieldHomogenization.Schemes.Proportional Type

Proportional() <: DifferentialTrajectory

All phases grow proportionally during the differential integration.

MeanFieldHomogenization.Schemes.Sequential Type

Sequential(order::Vector{Symbol}) <: DifferentialTrajectory

Introduce phases in the given order. The first phase ramps from 0 to its target fraction over the steps it owns, then is frozen; the next phase ramps over its allotted steps; and so on.

MeanFieldHomogenization.Schemes.CustomPath Type

CustomPath(path::Dict{Symbol, <:AbstractVector{<:Real}}) <: DifferentialTrajectory
CustomPath(:phase => values, ...)

Explicit per-phase trajectory. path[:phase] must be a length-N monotone vector with path[:phase][1] = 0 and path[:phase][end] = 1, where N is the number of differential steps.

The pair form is a convenience constructor for the same thing: CustomPath(:I1 => [0.0, 0.5, 1.0], :I2 => [0.0, 0.2, 1.0]).

MeanFieldHomogenization.Schemes.Path Type

Path(path::Dict{Symbol, <:Function}) <: DifferentialTrajectory
Path(:phase => τ -> f(τ), ...)

Explicit per-phase trajectory as a callable. path[:phase] is a function of the fictitious incorporation time τ ∈ [0, 1] returning the effective volume fraction ratio f_α(τ) / f_α^∞ ∈ [0, 1] for solid phases, or the density ratio ε_α(τ) / ε_α^∞ for crack phases — f(0) = 0, f(1) = 1, monotone non-decreasing.

The derivative df/dτ is computed by ForwardDiff.derivative at each ODE step. This is the natural API for the multi-phase DEM incorporation-sequence ODE :

with the volumetric balance dφ = (𝟙 - f ⊗ 𝐔)^{-1} · df (Sherman- Morrison) inverted at each τ to translate user-supplied f_α(τ) into the increments dφ_α / dτ.

The single-phase case is degenerate (f₁ itself serves as τ) and does not require a Path — the default Proportional is sufficient.

The pair form is a convenience constructor for the same thing:

Path(:I1 => τ -> τ^2, :I2 => τ -> 2τ - τ^2)

MeanFieldHomogenization.Schemes.AndersonDefault Type

AndersonDefault()

Marker selecting the built-in Anderson-accelerated fixed-point solver (default for SelfConsistent). Pure Julia, Dual-safe.

MeanFieldHomogenization.Schemes.NewtonDefault Type

NewtonDefault()

Marker selecting the built-in Newton-Raphson solver with ForwardDiff Jacobian (alternative to AndersonDefault).

MeanFieldHomogenization.Schemes.AutoNonlinear Type

AutoNonlinear()

Marker selecting an auto-resolving nonlinear solver: at each call it checks, at runtime, whether the weak extension MeanFieldHomogenizationNonlinearSolveExt is active (Base.get_extension) — i.e. whether NonlinearSolve.jl has been loaded into the session, whether by an explicit using NonlinearSolve or transitively through some other dependency (empirically, on a typical installation today the extension is already active as soon as MeanFieldHomogenization itself is loaded — the exact transitive path depends on the resolved dependency graph and is not guaranteed across versions).

• If active, it dispatches to a globalized SciML algorithm (NonlinearSolve.TrustRegion() — chosen over a plain NewtonRaphson() for its better robustness near the self-consistent bifurcation), through the same ForwardDiff-safe path (implicit-function-theorem lift) as any other NonlinearSolve.jl algorithm passed explicitly.

• Otherwise (a slimmed-down dependency set, or a future OrdinaryDiffEq that no longer needs NonlinearSolve.jl internally) it falls back to the dependency-free built-in NewtonDefault.

This gives a solver choice that always works, with or without NonlinearSolve.jl explicitly requested.

AutoNonlinear is not the default of SelfConsistent / AsymmetricSelfConsistent (which remains AndersonDefault): this is a numerical-robustness choice, not a dependency-cost one — the porous-percolation regime relies on the Picard positive-definite guard and select_best to track the physical branch through the bifurcation, a property a root-finder does not share. Pass algorithm = AutoNonlinear() explicitly to opt in.

Entry point

MeanFieldHomogenization.Core.homogenize Function

homogenize(cell, scheme, property::Symbol; kw...)

Compute the effective property of cell under scheme. Methods live in Schemes/homogenize.jl; the generic is declared here so that every cell type attaches to the same function.

MeanFieldHomogenization.Schemes.crack_family_compliances Function

crack_family_compliances(rve, scheme, C_hom; property = :C, kw...)

Decompose the macroscopic compliance of rve into a solid part and one contribution per crack family, such that

Returns a NamedTuple:

• solid — the solid part (for a single-matrix RVE, the matrix compliance );

• families — a Dict{Symbol,…} mapping each crack phase name to its , normalized per unit  , so that the scaled contribution is recovered with delta_compliance.

C_hom must be the estimate that scheme itself returned for this rve (homogenize(rve, scheme, property)): the decomposition is an identity of the scheme at its own fixed point, and evaluating it anywhere else is meaningless. crack_family_residual checks exactly that.

The are what drives an aperture update — see the derivation of   at the top of this file, and MeanFieldHomogenization.Constitutive, which consumes them.

Supported schemes

Scheme		correction
SelfConsistent
MoriTanaka		none (identity)

For Mori-Tanaka the identity    is exact — a two-line consequence of the 𝔹 : 𝔸⁻¹ body — but only when the matrix is the sole phase carrying volume. An RVE with solid inclusions and cracks is rejected rather than answered approximately.

Unsymmetrized families only

  is a statement about one orientation. A crack phase declared with symmetrize = :iso or :ti has already had its ℍ averaged over an orbit of orientations, and no single remains to project onto. Such a phase raises an ArgumentError here instead of returning a plausible but meaningless tensor. Discrete fracture families — the ARMA setting — are unsymmetrized by construction.

See also crack_family_residual, homogenize.

MeanFieldHomogenization.Schemes.crack_family_residual Function

crack_family_residual(rve, scheme, C_hom; property = :C, kw...) -> Real

Relative Frobenius residual of the identity guaranteed by crack_family_compliances:

A model self-check, and the regression test of the decomposition. For an iterative scheme the residual is bounded by the solver tolerance, not by machine epsilon: SelfConsistent returns step(x) on convergence, so with the default reltol = 1e-8 expect a residual of that order. Tighten abstol/reltol in the scheme to tighten the residual, and do not enable select_best (which returns an earlier iterate and so breaks the identity by that iterate's residual).

MeanFieldHomogenization.Schemes.fracture_permeability Function

fracture_permeability(k_matrix, cracks, densities; kw...) -> Tens{2,3}

Effective conductivity (permeability, diffusivity) of a solid of conductivity k_matrix crossed by families of flowing fractures.

• k_matrix — the matrix conductivity, a scalar or a Tens{2,3}. For a fractured rock it is small but must be non-zero: with an exactly impermeable matrix the medium conducts nothing below the percolation threshold, and the estimate degenerates.

• cracks — the ConductiveCrack families, each carrying its own normal and fracture conductivity.

• densities — the Budiansky crack density of each family.

Solves the self-consistent condition

each family being read in the effective medium, which is what lets the fractures see one another and produces a percolation threshold.

Keywords

abstol, reltol, maxiters, damping, verbose are passed to the fixed-point driver shared with SelfConsistent. Strong contrasts usually need damping < 1.

cracks = (ConductiveCrack(1.0; conductivity = 6.7e-11, euler_angles = (π/2, 0.0)),
          ConductiveCrack(1.0; conductivity = 6.7e-11, euler_angles = (π/2, π/2)))
K = fracture_permeability(1.0e-18, cracks, (0.37, 0.37))

Orientation lives in the crack, density in the argument

Same split as an RVE: the family's normal is part of its geometry, its abundance is not.

fracture_permeability(rve::RVE; solid = nothing, property = :K, kw...) -> Tens{2,3}

fracture_permeability reading the solid conductivity, the ConductiveCrack families and their densities straight off an RVE.

Phases that are not conductive cracks are ignored — a fracture network model has no use for them, and silently folding them in would be worse than saying so.

solid names the phase that plays the intact skeleton; left unset, it is the phase absorbing the volume complement.

MeanFieldHomogenization.Schemes.differential_path Function

differential_path(rve, scheme::DifferentialScheme, property::Symbol; kw...)
    -> (τ::Vector, states::Vector{<:AbstractTens})

Effective property of rve all along the differential scheme's fictitious incorporation time τ ∈ [0, 1], instead of only at τ = 1 as homogenize returns.

τ carries nsteps + 1 uniformly spaced save points (the saveat of the underlying ODE solve — the integration step itself remains adaptive and controlled by abstol / reltol), and states[k] is the effective property after the fraction of each phase has been grown to f_α(τ[k]) along the scheme's trajectory. states[1] is the matrix property and states[end] is exactly what homogenize returns.

Useful to plot the construction history of a composite, and to compare incorporation trajectories at equal target fractions.

τ, Cs = differential_path(rve, DifferentialScheme(; nsteps = 200), :C)
ks = [k_mu(C)[1] for C in Cs]

MeanFieldHomogenization.Schemes.SCHEME_ALIAS Constant

SCHEME_ALIAS :: Dict{Symbol, Type{<:HomogenizationScheme}}

Maps Symbol shortcuts to concrete scheme types for the convenience overload homogenize(rve, ::Symbol).

The canonical aliases are all lowercase (:voigt, :reuss, :dilute, :dilute_dual, :mori_tanaka, :mt, :maxwell, :pcw, :sc, :asc, :differential, :diff) for consistency with the algorithm-method symbols accepted elsewhere in the package (:auto, :residues, :decuhr, :nestedquadgk, :analytical). The CamelCase forms (:MoriTanaka, :Differential) and the ECHOES-compatible upper-case codes (:MT, :DIFF, …) are kept as extra aliases for ease of porting.

Symmetry projections

Best-fit projection of a tensor onto a symmetry class. These force major symmetry, unlike the exact rotation-group averages isotropify / transverse_isotropify (re-exported from TensND); the two differ whenever the input is not major-symmetric, and the difference is worked through in Symmetrization showcase.

best_fit_iso(t), best_fit_ti(t, axis) and best_fit_ortho(t, frame) live in TensND, next to the proj_tens machinery they wrap, and are re-exported here because turning a scheme result into publishable parameters is part of this package's surface. See the TensND documentation for their docstrings.

Orientation distributions

Discretizing a continuous orientation distribution into a finite set of phases, the alternative to declaring a symmetrize on a single one.

MeanFieldHomogenization.Schemes.polar_orientation_bins Function

polar_orientation_bins(N) -> Vector{@NamedTuple{θ, weight}}

Discretize the polar angle θ ∈ [0, π/2] into N families with solid-angle weights on the hemisphere, following the quasi-brittle strength model of [72] (2011) / echoes disc_theta convention :

• family angles θ_i = (π/2)·(i−1)/(N−1) (endpoints included),

• bin edges at the mid-points, θ_i^± = (π/2)·(i−1±1/2)/(N−1) clamped to [0, π/2],

• weights w_i = cos θ_i^− − cos θ_i^+ (Σ w_i = 1).

Declare each family as a phase with the tilted geometry (Spheroid(ω; euler_angles = (θ_i, 0, 0))), fraction f · w_i and symmetrize = TISymmetrize((0, 0, 1)) — the exact azimuthal average about the global axis represents the uniform-in-azimuth orbit of the family. As N → ∞ the family sum converges (in O(Δθ²)) to the single IsoSymmetrize phase.