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    over the other VolumeFraction values.

• 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 (Ponte Castañeda and Willis, 1995), 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 (Sevostianov and Giraud, 2013; Ponte Castañeda and Willis, 1995).

Future extension, not implemented: a PairwiseDistribution carrying a per-pair (i, j) ↦ shape mapping (Ponte Castañeda and Willis, 1995; Willis, 1982). 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 when the inclusions are spheres too. 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 concentration 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 concentration-tensor calculations without any RVE.

MeanFieldHomogenization.Schemes.NoSymmetrize Type

NoSymmetrize() <: AbstractSymmetrize

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

MeanFieldHomogenization.Schemes.IsoSymmetrize Type

IsoSymmetrize() <: AbstractSymmetrize

The concentration 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)) <: AbstractSymmetrize

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

The average is defined only when the reference medium of the auxiliary Eshelby problem is itself invariant under rotations about axis (isotropic, or transversely isotropic about axis): only then do averaging the concentration tensors and averaging the contributions give the same scheme. The reference is taken through the same average, which leaves such a reference unchanged; homogenize warns when it does change it.

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 a 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. fraction = :rest makes the phase absorb the volume complement — at most one phase per RVE, and not under closure = RescaledFractions().

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 concentration 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 embeds 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 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 concentration 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 as soon as the inclusions are spheres too (more generally, whenever every inclusion phase has the Hill tensor of the distribution shape). Supplying one silently therefore answered a scheme whose whole purpose is the 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 (     ). damping = 0.0 default; raise to about 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      (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 the weights of TensND's param_norms; 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,   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 close to ) 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   ,

2. Compute the Jacobian   via ForwardDiff.jacobian,

3. Take the Newton step    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 (at most 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) -> Vector

The weights of the Newton unknowns of the self-consistent solvers. The unknowns are the coordinates of TensND's to_vector (the coefficients of a structured class, the canonical Kelvin-Mandel components of a generic tensor), weighted by TensND's param_norms so that their Euclidean norm is the Frobenius norm of the tensor they encode: abstol and reltol then mean the same whatever the algorithm and the symmetry class. Newton is invariant under that diagonal rescaling (the root and the iterates are untouched); only the norm of the stopping test and of the Armijo condition changes, and a trust region is put in a physical metric.

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 , 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 the flat-inclusion Hill tensor of layer and     . Serves elasticity and transport, and accepts imperfect interfaces of spring / membrane / Kapitza / surface-conductive type. See the theory page on the laminate and Backus (1962).

MeanFieldHomogenization.Schemes.Dilute Type

Dilute() <: HomogenizationScheme

Dilute scheme:     where     is the size-independent stiffness contribution (Eshelby, 1957; Kachanov and Sevostianov, 2018).

MeanFieldHomogenization.Schemes.DiluteDual Type

DiluteDual() <: HomogenizationScheme

Dual dilute scheme on the compliance:    , returning  .

MeanFieldHomogenization.Schemes.MoriTanaka Type

MoriTanaka() <: HomogenizationScheme

Mori-Tanaka scheme (Mori and Tanaka, 1973; Christensen, 1990).

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   and a stiff one at  , a Voigt seed and a phase seed agree to 1e-9 or better. They part only next to the porous percolation threshold ( , ), 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 Jacobian is exactly singular. Hence the default above.

The seed also sets a floor

Every solver keeps each component of its iterate above times the largest component of the seed (and above ): the positive-definiteness guard that keeps the iteration off the null fixed point near percolation. A solution with a smaller component is out of reach — the iteration stops on the floor, and the result is wrong without notice. With the default seed of a solid around 1, that floor is about 1e-8. When the moduli of the RVE lie far below the seed — soft interfaces, a loose granular assembly — seed with a tensor of the order of the expected solution, for instance init = iso_stiffness(kt, kt).

Standard kwargs forwarded to the solver: abstol, reltol, maxiters, damping, verbose, select_best. Convergence is declared when     (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     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   (Norris, 1985) :

with the volume balance     inverted by Sherman- Morrison so the user supplies effective volume fractions 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 the SciML solve).

• alg — explicit ODE algorithm. nothing selects Tsit5() (5th order adaptive Runge-Kutta), the only one MeanFieldHomogenization loads. Any other algorithm of the OrdinaryDiffEq family is accepted once its package is loaded in the session: using OrdinaryDiffEq for all of them, or one family alone (using OrdinaryDiffEqVerner for Vern9(), using OrdinaryDiffEqRosenbrock for Rodas5P()). Implicit algorithms (Rosenbrock23(), Rodas5P()) work with their default ForwardDiff Jacobian and time gradient: the right-hand side accepts dual numbers in the state and in τ.

• kwargs... — any other keyword is forwarded verbatim to the SciML 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  , 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  .

• Path — explicit per-phase trajectory as a callable τ -> f(τ) (auto-differentiated by ForwardDiff); the natural API for the multi-phase incorporation-sequence ODE (Norris, 1985).

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   returning the effective volume fraction ratio   for solid phases, or the density ratio for crack phases —  ,  , monotone non-decreasing.

The derivative 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     (Sherman- Morrison) inverted at each to translate user-supplied into the increments .

The single-phase case is degenerate ( 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, by an explicit using NonlinearSolve or through another package that loads it. MeanFieldHomogenization itself does not: since 0.20.2 it depends on OrdinaryDiffEqTsit5 alone, which does not load NonlinearSolve.jl (the full OrdinaryDiffEq did, which made the extension active as soon as the package was loaded).

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

with the Budiansky density of crack family .

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

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  , instead of only at   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.

Concentration tensors of an estimate and Levin's theorem

MeanFieldHomogenization.Schemes.ConcentrationTensors Type

ConcentrationTensors

Concentration tensors of every phase of a cell under one scheme, as returned by concentration_tensors.

For each phase r, the strain concentration tensor and the stress–strain one map the macroscopic strain onto the mean strain and the mean stress of the phase,    and   . They satisfy

being the volume fraction of the phase. A crack family has no volume: its weight is its density, its strain tensor is the strain carried by its openings per unit density, and its stress tensor is zero. In conduction the same holds with the gradient for the strain and the opposite of the flux for the stress; the fields keep the names strain and stress for both orders.

Read the tensors with gradient_gradient_concentration, flux_gradient_concentration, gradient_flux_concentration and flux_flux_concentration, the contributions with property_contribution and inverse_property_contribution, and apply Levin's theorem with levin_free_gradient and levin_preflux: the same functions in elasticity and in conduction, whose names of elasticity (strain_strain_concentration, stiffness_contribution, levin_free_strain, …) are aliases for fourth-order estimates.

Fields

• effective — (or );

• compliance — its inverse;

• names — the phases, in declaration order;

• strain, stress — the tensors of each phase;

• weights — the volume fraction of each phase, or the density of a crack family;

• cracks — the crack families among the phases;

• properties — the property of each homogeneous phase, usable as a reference;

• reference — the default reference phase of the contributions, nothing when the scheme has none.

MeanFieldHomogenization.Schemes.concentration_tensors Function

concentration_tensors(cell, scheme, property::Symbol; kw...) -> ConcentrationTensors
concentration_tensors(cell, scheme; property = :C, kw...)

Concentration tensors of every phase of cell under scheme, for the property property (:C in elasticity, :K in conduction), together with the effective tensor they are consistent with. The keyword arguments are those of homogenize.

Every scheme on an RVE is supported. For a CrackDensity phase the tensors are given per unit density; the bounds ignore cracks, as they do in homogenize.

ct = concentration_tensors(rve, MoriTanaka(:M), :C)
A = gradient_gradient_concentration(ct, :I)        # ⟨ε⟩_I = A : E
N = property_contribution(ct, :I)             # ⟨σ − C_M : ε⟩_I = N : E

MeanFieldHomogenization.Schemes.levin_free_gradient Function

levin_free_gradient(ct::ConcentrationTensors, eigengradients) -> AbstractTens

Macroscopic free gradient of a composite whose phases carry uniform eigengradients, by Levin's theorem.

In elasticity the eigengradients are eigenstrains (free strains) : with     in each phase,     with

being the flux (stress) concentration tensor of the phase (flux_flux_concentration). In conduction they are vectors , the gradients at which the phases conduct no flux,     , and the composite obeys      with    .

eigengradients maps phase names to second-order tensors in elasticity, to vectors (first-order tensors) in conduction; phases left out carry none. A crack family holds no material, so an eigengradient in it has no effect and is refused. Explicit name in elasticity: levin_free_strain.

MeanFieldHomogenization.Schemes.levin_preflux Function

levin_preflux(ct::ConcentrationTensors, prefluxes) -> AbstractTens

Macroscopic pre-flux of a composite whose phases carry uniform pre-fluxes, the fluxes they carry under a vanishing gradient, by Levin's theorem.

In elasticity the pre-fluxes are prestresses : with     in each phase,     with

being the gradient (strain) concentration tensor of the phase. In conduction,      gives      with    .

A crack family enters with its density as weight and its tensor per unit density, and only the load on its faces matters. A pore pressure is the prestress   of the pores and of the cracks it fills, and the Biot tensor follows as       . Explicit name in elasticity: levin_prestress.

MeanFieldHomogenization.Schemes.levin_free_strain Function

levin_free_strain(ct::ConcentrationTensors, eigenstrains) -> AbstractTens{2,3}

Name in elasticity of levin_free_gradient: the macroscopic free strain from the eigenstrains of the phases, for a fourth-order estimate only.

MeanFieldHomogenization.Schemes.levin_prestress Function

levin_prestress(ct::ConcentrationTensors, prestresses) -> AbstractTens{2,3}

Name in elasticity of levin_preflux: the macroscopic prestress from the prestresses of the phases, for a fourth-order estimate only.

The readers of a ConcentrationTensors are methods of gradient_gradient_concentration, flux_gradient_concentration, gradient_flux_concentration, flux_flux_concentration, property_contribution and inverse_property_contribution, for both orders, and of their names in elasticity and in conduction, documented in API — Concentration.

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 project 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   into families with solid-angle weights on the hemisphere, following the quasi-brittle strength model of Pichler and Hellmich (2011) and the disc_theta convention of Echoes:

• family angles    (endpoints included),

• bin edges at the mid-points,    clamped to ,

• weights    ( ).

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