API — Viscoelasticity

MeanFieldHomogenization.Viscoelasticity Module

MeanFieldHomogenization.Viscoelasticity

Ageing linear viscoelastic (ALV) homogenization. Provides:

• ViscoLaw — relaxation R(t,t') or creep J(t,t') kernel, scalar- or 4-tensor-valued, with built-in Maxwell / Kelvin constructors.

• trapezoidal_matrix — discretization of the Stieltjes integral on a time grid into a lower-block-triangular matrix (n×n for scalar, 6n×6n for 4-tensor in Mandel form).

• volterra_inverse — block forward-substitution that takes a discrete relaxation kernel to the corresponding creep kernel (and vice versa).

• iso_params_from_blocks / iso_blocks_from_params and their ti_ / ortho_ counterparts — conversions between symmetry-structured per-component scalar matrices and the full 6n×6n block matrix.

• hill_kernel — discrete ALV Hill polarization tensor for an ellipsoidal inclusion, isotropic-matrix branch using the time-space decoupling formula ([9], App. ALV Hill kernel).

• Time-domain viscoelastic homogenization schemes (Voigt, Reuss, Dilute, DiluteDual, Mori-Tanaka, Maxwell, Self-Consistent), plugged into the existing homogenize dispatcher whenever a phase carries a ViscoLaw property.

All ALV operators are stored as dense Matrix{T} of size (B·n)×(B·n) (B = 6 for 4-tensor, B = 1 for scalar) with explicit zeros above the block diagonal — this is the convention of [60] and the C++ ECHOES reference.

Constitutive laws

MeanFieldHomogenization.Viscoelasticity.AbstractViscoLaw Type

AbstractViscoLaw

Root supertype for viscoelastic kernels. Concrete subtypes provide an eval(law, t, t_p) method returning a scalar or a 4-tensor.

MeanFieldHomogenization.Viscoelasticity.ViscoLaw Type

ViscoLaw(eval_fun, mode::Symbol = :relaxation)

Concrete viscoelastic law. Wraps an eval_fun::F callable (t, t_p) -> X together with its mode Symbol (:relaxation or :creep). The output type X may be:

• a Real or Complex scalar (scalar Volterra kernel) ;

• a TensND.AbstractTens{4,3} (4-tensor relaxation / creep tensor) ;

• an AbstractMatrix of size 6×6 already in Mandel form.

Construct via ViscoLaw directly or use the convenience constructors maxwell_relaxation, kelvin_creep, maxwell_iso, kelvin_iso.

MeanFieldHomogenization.Viscoelasticity.visco_mode Function

visco_mode(law) -> Symbol

Return the mode (:relaxation or :creep) of law.

MeanFieldHomogenization.Viscoelasticity.VALID_VISCO_MODES Constant

VALID_VISCO_MODES :: Tuple{Symbol, Symbol}

Allowed values for ViscoLaw.mode :

• :relaxation — R(t, t') maps strain history to stress.

• :creep — J(t, t') maps stress history to strain.

MeanFieldHomogenization.Viscoelasticity.visco_eval Function

visco_eval(law, t, t_p)

Evaluate the kernel at (t, t_p).

MeanFieldHomogenization.Viscoelasticity.maxwell_relaxation Function

maxwell_relaxation(C_inf, C_branches, taus; mode = :relaxation)

Build a generalized Maxwell relaxation kernel

R(t, t') = C_inf + Σ_i C_branches[i] · exp(-(t - t')/taus[i])

for t ≥ t', 0 otherwise. C_inf and C_branches[i] may be scalars or 4-tensors (TensND.AbstractTens{4,3}); taus[i] is a positive relaxation time of the i-th branch.

MeanFieldHomogenization.Viscoelasticity.kelvin_creep Function

kelvin_creep(J_0, J_branches, taus; mode = :creep)

Build a Kelvin (or Kelvin-Voigt-Generalized) creep kernel

J(t, t') = J_0 + Σ_i J_branches[i] · (1 - exp(-(t - t')/taus[i]))

for t ≥ t', 0 otherwise. J_0 is the instantaneous compliance, J_branches[i] the i-th branch compliance, taus[i] its retardation time.

MeanFieldHomogenization.Viscoelasticity.maxwell_iso Function

maxwell_iso(k, mu, eta_k, eta_mu) -> ViscoLaw

Convenience: build an isotropic Maxwell relaxation 4-tensor kernel

R(t, t') = 3 k · exp(-(t - t')/eta_k) · 𝕁
         + 2 mu · exp(-(t - t')/eta_mu) · 𝕂

with 𝕁 the spherical projector (1/3) 𝟙 ⊗ 𝟙 and 𝕂 = 𝕀 - 𝕁 the deviatoric projector. The output is a TensND.TensISO{4,3} at every (t, t') with t ≥ t'.

MeanFieldHomogenization.Viscoelasticity.kelvin_iso Function

kelvin_iso(k_0, mu_0, k_branches, mu_branches, taus_k, taus_mu) -> ViscoLaw

Convenience: build an isotropic Kelvin creep 4-tensor kernel

J(t, t') = (1/(3 k_0)) 𝕁 + (1/(2 mu_0)) 𝕂
         + Σ_i (1/(3 k_branches[i])) (1 - exp(-(t-t')/taus_k[i])) 𝕁
         + Σ_i (1/(2 mu_branches[i])) (1 - exp(-(t-t')/taus_mu[i])) 𝕂

k_branches, mu_branches, taus_k, taus_mu may be empty if no Kelvin branches are required (instantaneous-only compliance).

MeanFieldHomogenization.Viscoelasticity.heaviside_law Function

heaviside_law(C; mode = :relaxation)

Return a ViscoLaw corresponding to a purely elastic kernel R(t,t') = C · H(t-t'). Useful for testing the elastic limit of the ALV pipeline.

Trapezoidal discretization

MeanFieldHomogenization.Viscoelasticity.trapezoidal_matrix Function

trapezoidal_matrix(law::ViscoLaw, times::AbstractVector{<:Real}) -> Matrix

Build the discrete (B·n) × (B·n) lower-block-triangular matrix representing the Stieltjes integral    on the time grid times = (t_0, …, t_{n-1}) using the trapezoidal rule of [60]. B = 1 for scalar-valued kernels and B = 6 for 4-tensor- or 6×6-matrix-valued kernels (Mandel form).

The exact block layout is described at the top of trapezoidal.jl. The output is type-stable in the element type returned by visco_eval(law, ...).

MeanFieldHomogenization.Viscoelasticity._trapezoidal_relaxation Function

_trapezoidal_relaxation(law::ViscoLaw, times, B) -> Matrix

Build the discrete relaxation block matrix from a ViscoLaw, regardless of whether the law is in :relaxation or :creep mode. When the law is :creep, the trapezoidal compliance matrix is inverted (block forward substitution at block_size = B) to obtain the corresponding relaxation matrix — the convention every ALV scheme assumes internally.

B is the block size (6 for order-4 4-tensor / Mandel, 3 for order-2 vector-tensor).

Volterra algebra

MeanFieldHomogenization.Viscoelasticity.volterra_inverse Function

volterra_inverse(M::AbstractMatrix; block_size::Int = 6) -> Matrix

Compute the Volterra inverse of a lower-block-triangular matrix M of size (B·n) × (B·n) with B = block_size. The result is also lower-block-triangular and satisfies M * volterra_inverse(M) = H 𝟙 (block-diagonal identity).

block_size must be a positive divisor of size(M, 1) and size(M, 2), typically 1 (scalar Volterra kernel) or 6 (4-tensor in Mandel form).

The cost is O(B³ n²) flops via block forward-substitution; each diagonal block is inverted with the dense inv(...) of stdlib (so B = 6 costs only a few hundred flops per block).

volterra_inverse(K::AbstractALVKernel) -> AbstractALVKernel

Volterra inverse of a structured ALV kernel. Stays in the same symmetry class : iso ↦ iso, TI ↦ TI, ortho ↦ ortho. Avoids materializing the (6n × 6n) matrix.

MeanFieldHomogenization.Viscoelasticity.volterra_product Function

volterra_product(A::AbstractMatrix, B::AbstractMatrix) -> Matrix

Discrete Volterra product A ∘ B: a regular matrix multiplication of two lower-block-triangular matrices, returned as a fresh Matrix. The result is again lower-block-triangular by construction.

Equivalent to A * B, exposed as a named function to make the viscoelastic algebra explicit at call sites.

MeanFieldHomogenization.Viscoelasticity.volterra_divide Function

volterra_divide(M, S; block_size = 1) -> Matrix

Compute T = M ∘ S^{-vol} (Volterra-divide) by direct block forward substitution on the linear system T · S = M. Numerically stable where M * volterra_inverse(S) would lose precision through a huge intermediate inverse, e.g. for soft-phase moduli (κ, μ → 0) or step-activated ViscoLaw kernels in multi-layer ALV recurrences.

block_size must divide size(M, 1); typical values are 1 (scalar Volterra) and 6 (4-tensor Mandel). The result is lower-block-triangular if both M and S are.

MeanFieldHomogenization.Viscoelasticity.volterra_left_divide Function

volterra_left_divide(S, M; block_size = 1) -> Matrix

Compute T = S^{-vol} ∘ M (Volterra LEFT-divide) by direct forward substitution on the linear system S · T = M.

Use this rather than volterra_divide when the closed-form algebra requires the inverse on the LEFT of the numerator (e.g. the Hervé–Zaoui bulk/shear transition matrices, where T = M_b^{-1} · M_a). On a non-uniform time grid Volterra trapezoidal matrices do not form a commutative algebra, so the order matters and right-vs-left divides give different results.

block_size must divide size(M, 1); typical values are 1 (scalar Volterra) and 6 (4-tensor Mandel). Both arguments must be lower- block-triangular of the same size; the result is lower-block-triangular.

volterra_left_divide(S::AbstractALVKernel, M::AbstractALVKernel)
    -> AbstractALVKernel

Volterra left-divide T = S^{-vol} ∘ M within the structured class. Auto-promotes mixed inputs (e.g. iso S + TI M) to the more general class before solving.

Symmetry-class conversions

MeanFieldHomogenization.Viscoelasticity.iso_params_from_blocks Function

iso_params_from_blocks(M) -> (α::Matrix, β::Matrix)

Decompose a 6n×6n block matrix whose every 6×6 block is an isotropic 4-tensor in Mandel form into the two scalar parameter matrices α and β, both of size n × n. The block layout is

  M_block(i,j) = α[i,j] · 𝕁_Mandel + β[i,j] · 𝕂_Mandel

where 𝕁_Mandel = (1/3) e₁ e₁ᵀ (e₁ the Mandel-1 unit) augmented with the symmetric off-diagonal part (1/3) for the upper 3×3 block, and 𝕂_Mandel = 𝕀_Mandel - 𝕁_Mandel.

α[i,j] and β[i,j] are extracted from the diagonal entry M_block[1,1] and the (4,4) Mandel-shear entry of each block: α[i,j] = (M_block[1,1] + 2 M_block[1,2]) β[i,j] = M_block[4,4].

MeanFieldHomogenization.Viscoelasticity.iso_blocks_from_params Function

iso_blocks_from_params(α::Matrix, β::Matrix) -> Matrix

Inverse of iso_params_from_blocks: build a 6n×6n block matrix whose every block is the iso 4-tensor α[i,j] · 𝕁 + β[i,j] · 𝕂 in Mandel form.

Both α and β must be n × n and have a common element type.

Note: a kron(α, 𝕁_M) + kron(β, 𝕂_M) formulation was experimented (P3.3.1) and reverted — the two extra (6n × 6n) intermediate allocations made it ~3× slower than this single-pass scalar loop, which writes each entry exactly once.

MeanFieldHomogenization.Viscoelasticity.ti_params_from_blocks Function

ti_params_from_blocks(M; axis = (0, 0, 1))
    -> NTuple{6, Matrix{T}}

Decompose a 6n×6n block matrix whose every 6×6 block is a TI 4-tensor with axis n (only n = e₃ is currently supported) into the six scalar Walpole parameter matrices (ℓ₁, ℓ₂, ℓ₃, ℓ₄, ℓ₅, ℓ₆), each of size n × n.

MeanFieldHomogenization.Viscoelasticity.ti_blocks_from_params Function

ti_blocks_from_params(ℓ::NTuple{6, AbstractMatrix}; axis = (0, 0, 1))
    -> Matrix{T}

Inverse of ti_params_from_blocks: rebuild a 6n×6n block matrix whose every 6×6 block has the TI Mandel structure with axis n = e₃ and Walpole coefficients (ℓ₁[i,j], …, ℓ₆[i,j]).

MeanFieldHomogenization.Viscoelasticity.ortho_params_from_blocks Function

ortho_params_from_blocks(M; axes = ((1,0,0),(0,1,0),(0,0,1)))
    -> NTuple{12, Matrix{T}}

Decompose a 6n × 6n block matrix whose every 6×6 block is an ortho 4-tensor with material frame (e₁, e₂, e₃) (only the canonical frame is currently supported) into the 12 scalar parameter matrices

(o₁₁, o₁₂, o₁₃, o₂₁, o₂₂, o₂₃, o₃₁, o₃₂, o₃₃, o₄, o₅, o₆),

each of size n × n, where the 3×3 normal block of the Mandel form is stored row-major in entries 1..9 and the 3 shear-diagonal Mandel entries (M[4,4], M[5,5], M[6,6]) map to (o₄, o₅, o₆).

Note that closure under Volterra product does not preserve major symmetry (o₁₂ ≠ o₂₁ in general), so the full 9-entry normal block is needed even for major-symmetric inputs.

MeanFieldHomogenization.Viscoelasticity.ortho_blocks_from_params Function

ortho_blocks_from_params(o::NTuple{12, AbstractMatrix};
                          axes = ((1,0,0),(0,1,0),(0,0,1)))
    -> Matrix{T}

Inverse of ortho_params_from_blocks: rebuild a 6n × 6n block matrix whose every 6×6 block has the ortho Mandel structure with axes (e₁, e₂, e₃). All off-block-diagonal entries are zero (block 1..3 ↔ block 4..6 coupling is forbidden by orthotropic symmetry).

Hill ALV kernel

MeanFieldHomogenization.Viscoelasticity.hill_kernel Function

hill_kernel(ell::AbstractEllipsoidalInclusion,
            C0_law::ViscoLaw,
            times::AbstractVector{<:Real}) -> Matrix

Discrete ALV Hill polarization tensor kernel for the inclusion ell in an isotropic ALV matrix described by C0_law (a ViscoLaw returning a TensND.TensISO{4,3} 4-tensor at each (t, t')).

The output is a (6n × 6n) lower-block-triangular Matrix{T} with n = length(times), in the same Mandel block convention as trapezoidal_matrix.

Implementation follows the time-space decoupling formula of the ECHOES manual appendix viscoelastic_hill_kernel.qmd:

1. Discretize the matrix kernel C0_law on times to a 6n × 6n block matrix R̃_M.

2. Extract the iso scalar parameter matrices (α, β) of R̃_M (α = 3K, β = 2μ per Mandel convention).

3. Build the longitudinal M_long = (α + 2β)/3 and the shear M_shear = β/2 n×n matrices.

4. Take the scalar Volterra inverses J_long = M_long^{-vol} and J_shear = M_shear^{-vol}.

5. Compute the Mandel forms U^A_M and V^A_M of the elastic auxiliary tensors tens_UA(ell) and tens_VA(ell).

6. Assemble block-by-block: P̃_block(i,j) = J_long[i,j] · U^A_M + J_shear[i,j] · (V^A_M - U^A_M).

Schemes — generic (6n × 6n) algebra

MeanFieldHomogenization.Viscoelasticity.voigt_alv Function

voigt_alv(C_phases::AbstractVector, fractions::AbstractVector) -> Matrix

Voigt (uniform-strain) bound: C_eff = Σ_r f_r · C̃^r. Each C_phases[r] is a (6n × 6n) block matrix and fractions[r] is the volume fraction of phase r.

MeanFieldHomogenization.Viscoelasticity.reuss_alv Function

reuss_alv(C_phases::AbstractVector, fractions::AbstractVector) -> Matrix

Reuss (uniform-stress) bound: invert each C̃^r, do the volume average of compliances, then invert the result.

J̃_eff = Σ_r f_r · J̃^r , C̃_eff = (J̃_eff)^

MeanFieldHomogenization.Viscoelasticity.dilute_alv Function

dilute_alv(C_0, contribs, fractions) -> Matrix

Dilute scheme: C̃_eff = C̃^0 + Σ_r f_r · Ñ^{r,dil} where the dilute contributions Ñ^{r,dil} of each inclusion phase are pre-computed via dilute_contribution_alv (each (6n × 6n)), f_r is the volume fraction of phase r, and the matrix C̃^0 is the reference.

Note that for the dilute scheme the matrix volume fraction f_0 does not appear: the inclusions are treated as if they were embedded in an infinite matrix.

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv Function

dilute_dual_alv(C_0, contribs, fractions) -> Matrix

DiluteDual scheme: the dilute contribution is computed in compliance space. Equivalent in form to dilute_alv but with the matrix acting on the compliance kernel. Useful when the inclusions are weak (e.g. cracks) — the linearization around the matrix compliance converges better.

J̃_eff = J̃^0 + Σ_r f_r · H̃^{r,dil} C̃_eff = (J̃_eff)^

with H̃^{r,dil} = (J^r - J^0) ∘ B̃^{r,dil} and B̃^{r,dil} = (𝟙 + (𝟙 - Ã^{r,dil}))^{-vol} … see Sevostianov 2008 for the dual formulation in the elastic case.

Implementation routes through volterra_inverse of C_0 to obtain J^0 and back.

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv Function

mori_tanaka_alv(C_0, A_duts, contribs, fractions, f_matrix) -> Matrix

Mori-Tanaka scheme: the average matrix-strain is taken as the reference "infinite-matrix" strain, hence each inclusion experiences a perturbed strain proportional to the inverse of the volume-weighted concentration average.

C̃_eff = C̃^0 + (Σ_r f_r Ñ^{r,dil}) ∘ (f_0 · 𝟙 + Σ_s f_s Ã^{s,dil})^

A_duts[r] is the dilute concentration kernel of phase r, contribs[r] = ΔC^r ∘ A_duts[r], fractions[r] its volume fraction in the RVE, and f_matrix the matrix volume fraction f_0.

MeanFieldHomogenization.Viscoelasticity.maxwell_alv Function

maxwell_alv(C_0, contribs, fractions; H_0) -> Matrix

Maxwell scheme: the reference is replaced by a "host medium" with a prescribed distribution shape, whose Hill kernel is H_0 = P̃_d. The Volterra-discrete formula is

C̃_eff = C̃^0 + Σ̃ ∘ (𝟙 - P̃_d ∘ Σ̃)^

where Σ̃ = Σ_r f_r Ñ^{r,dil} and the dilute contributions are the ones computed by dilute_contribution_alv.

When the distribution shape coincides with the inclusion shape and the matrix is the reference, Maxwell reduces to Mori-Tanaka.

Pass H_0 = hill_kernel(distribution_shape, C_0_law, times) for the default uniform-sphere distribution.

MeanFieldHomogenization.Viscoelasticity.self_consistent_alv Function

self_consistent_alv(rve, prop; times,
                    abstol = 1e-10, reltol = 1e-8, maxiters = 200,
                    damping = 0.0, verbose = false,
                    select_best = false) -> Matrix

Self-consistent ALV homogenization. Iterates the symmetric Picard fixed point on the (6n × 6n) block matrix until convergence.

The initial estimate is the discretized matrix kernel C̃^0. Each iteration rebuilds the per-phase Hill kernels using the current estimate's iso parameters, computes the dilute concentration tensors, and forms C̃_{m+1}.

Returns the converged effective relaxation matrix.

Keyword arguments

• abstol — absolute Frobenius tolerance on ‖C̃_{m+1} − C̃_m‖.

• reltol — additive relative tolerance (multiplied by ‖C̃_m‖).

• maxiters — hard iteration cap.

• damping — Picard relaxation 0 ≤ damping < 1 (0 = no damping).

• verbose — print residual norms each iteration.

• select_best— return the best iterate seen (rather than the last) when convergence stalls.

MeanFieldHomogenization.Viscoelasticity.asymmetric_self_consistent_alv Function

asymmetric_self_consistent_alv(rve::RVE, prop::Symbol; times,
                                abstol = 1e-10, reltol = 1e-8,
                                maxiters = 200, damping = 0.0,
                                verbose = false, select_best = false)
    -> Matrix{T}

Asymmetric self-consistent viscoelastic homogenization. The iteration update reads

`C^{n+1} = C_M + Σ_i f_i (C_i − C_M) ∘ A^{dil,i}(C^n)`,

mirroring the reference ASC form. Returns the (6n × 6n) effective relaxation matrix once the residual ‖C^{n+1} − C^n‖_F falls below abstol + reltol · ‖C^n‖_F (or after maxiters Picard steps).

MeanFieldHomogenization.Viscoelasticity.pcw_alv Function

pcw_alv(C_0, contribs, fractions; H_dist) -> Matrix

Ponte-Castañeda & Willis (1995) viscoelastic homogenization in single-distribution-shape form. The formula is algebraically identical to maxwell_alv ; the only difference is that the Hill kernel H_dist is computed against the distribution shape stored in rve.distribution_shape, not against any individual phase.

Use homogenize_alv with scheme = PonteCastanedaWillis() — the dispatcher reads the RVE's distribution shape and forwards here.

MeanFieldHomogenization.Viscoelasticity.differential_alv Function

differential_alv(rve::RVE, prop::Symbol; times,
                  nsteps = 100, trajectory = nothing,
                  abstol = 1e-8, reltol = 1e-6, alg = nothing,
                  formulation = :stiffness) -> Matrix{T}

Differential homogenization in ageing linear viscoelasticity, solved as a SciML ODE on the fictitious incorporation time τ ∈ [0, 1] ([31]; user's hand-written DEM note) :

with the volume balance df = (𝟙 − f ⊗ 𝐔)·dφ inverted by Sherman- Morrison for solid phases (cracks contribute their density derivative directly). nsteps is the density of save points along τ ; the integration step is controlled by abstol / reltol.

formulation = :compliance integrates the dual ODE on the creep function J̃ = C̃^{-vol} instead, through H̃_α = −J̃ ∘ Ñ_α ∘ J̃, and inverts the result — the same choice as the elastic DifferentialScheme.

Supported inclusion geometries: ellipsoids / spheroids (through the ALV Hill kernel), LayeredSphere (through the bulk + shear ALV recurrences) and flat cracks (CrackDensity). Because the reference of the ODE is the running medium and the ALV Hill kernel is isotropic only, every phase must keep that medium isotropic — see _alv_diff_keeps_iso; the ODE throws an explicit ArgumentError otherwise instead of returning a wrong answer.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv Function

dilute_concentration_alv(C_E, C_0, P) -> Matrix

Dilute strain concentration kernel Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol} ([59] eq. 16). All inputs are (6n × 6n) block matrices ; the result is also (6n × 6n) and lower-block-triangular.

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv Function

dilute_contribution_alv(C_E, C_0, P) -> Matrix

Dilute strain contribution kernel Ñ = ΔC̃ ∘ Ã^dil ([59] eq. 17). This is the size-independent stiffness contribution of a single inclusion.

Internal guards

MeanFieldHomogenization.Viscoelasticity._alv_diff_keeps_iso Function

_alv_diff_keeps_iso(geom, sym, C_r) -> Bool

Whether a solid phase leaves the running effective medium of the differential ALV ODE inside the isotropic class.

This matters only for the differential scheme. The ALV Hill kernel hill_kernel exists for an isotropic reference only, and where Mori-Tanaka or the dilute scheme evaluate it against the (fixed, isotropic) matrix, the differential scheme evaluates it against the running estimate C̃(τ) — which an aligned, non-spherical inclusion progressively takes out of the iso class. Rather than silently reading (α, β) off a matrix that is no longer isotropic, the ODE refuses to start (see differential_alv).

An isotropic orientation average (symmetrize = :iso) restores the property for any shape, which is what randomly oriented inclusions or cracks mean physically.

Schemes — iso fast path

MeanFieldHomogenization.Viscoelasticity.voigt_alv_iso Function

voigt_alv_iso(αβ_phases, fractions) -> (α_eff, β_eff)

Iso-form Voigt bound: α_eff = Σ_r f_r α_r, β_eff = Σ_r f_r β_r.

MeanFieldHomogenization.Viscoelasticity.reuss_alv_iso Function

reuss_alv_iso(αβ_phases, fractions) -> (α_eff, β_eff)

Iso-form Reuss bound: invert per-phase compliances, average, invert back — done independently on the two scalar Volterra components.

MeanFieldHomogenization.Viscoelasticity.dilute_alv_iso Function

dilute_alv_iso(αβ_0, contribs_iso, fractions) -> (α_eff, β_eff)

Iso-form Dilute scheme: αβ_eff = αβ_0 + Σ_r f_r · (Ñ_r in iso form).

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_iso Function

dilute_dual_alv_iso(αβ_0, contribs_compliance_iso, fractions) -> (α_eff, β_eff)

Iso-form DiluteDual: invert the matrix to compliance space, average, invert back.

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_iso Function

mori_tanaka_alv_iso(αβ_0, A_duts_iso, contribs_iso, fractions, f_M)
    -> (α_eff, β_eff)

Iso-form Mori-Tanaka: C̃_eff = C̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}, reduced to two scalar Volterra problems on (α, β).

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_iso Function

maxwell_alv_iso(αβ_0, contribs_iso, fractions, αβ_H_0) -> (α_eff, β_eff)

Iso-form Maxwell scheme: C̃_eff = C̃_0 + Σ̃ ∘ (𝟙 - P̃_d ∘ Σ̃)^{-vol}, reduced to two scalar Volterra problems on (α, β).

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_iso Function

dilute_concentration_alv_iso(αβ_E, αβ_0, αβ_P) -> (α_A, β_A)

Iso-form dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol} reduced to two scalar Volterra problems on (α_E − α_0, β_E − β_0).

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_iso Function

dilute_contribution_alv_iso(αβ_E, αβ_0, αβ_P) -> (α_N, β_N)

Iso-form dilute contribution Ñ = ΔC̃ ∘ Ã^dil.

Schemes — TI fast path

MeanFieldHomogenization.Viscoelasticity.voigt_alv_ti Function

voigt_alv_ti(ℓ_phases, fractions) -> NTuple{6, Matrix}

TI-form Voigt bound: ℓᵢ_eff = Σ_r f_r ℓᵢ_r for each Walpole component.

MeanFieldHomogenization.Viscoelasticity.reuss_alv_ti Function

reuss_alv_ti(ℓ_phases, fractions) -> NTuple{6, Matrix}

TI-form Reuss bound: invert per-phase compliances in TI Walpole form, average, invert back.

MeanFieldHomogenization.Viscoelasticity.dilute_alv_ti Function

dilute_alv_ti(ℓ_0, contribs_ti, fractions) -> NTuple{6, Matrix}

TI-form Dilute scheme: ℓ_eff = ℓ_0 + Σ_r f_r · Ñ_r (component-wise).

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_ti Function

dilute_dual_alv_ti(ℓ_0, contribs_compliance_ti, fractions)
    -> NTuple{6, Matrix}

TI-form DiluteDual: invert to compliance Walpole form, average, invert back.

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_ti Function

mori_tanaka_alv_ti(ℓ_0, A_duts_ti, contribs_ti, fractions, f_M)
    -> NTuple{6, Matrix}

TI-form Mori-Tanaka: C̃_eff = C̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}, all in TI Walpole form.

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_ti Function

maxwell_alv_ti(ℓ_0, contribs_ti, fractions, ℓ_H_0) -> NTuple{6, Matrix}

TI-form Maxwell scheme.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_ti Function

dilute_concentration_alv_ti(ℓ_E, ℓ_0, ℓ_P) -> NTuple{6, Matrix}

TI-form dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol} reduced to a (2n)×(2n) block-Volterra inverse + two n×n scalar inverses.

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_ti Function

dilute_contribution_alv_ti(ℓ_E, ℓ_0, ℓ_P) -> NTuple{6, Matrix}

TI-form dilute contribution Ñ = ΔC̃ ∘ Ã^dil.

Schemes — ortho fast path

MeanFieldHomogenization.Viscoelasticity.voigt_alv_ortho Function

voigt_alv_ortho(o_phases, fractions) -> NTuple{12, Matrix}

Ortho-form Voigt bound: oₖ_eff = Σ_r f_r oₖ_r for each component.

MeanFieldHomogenization.Viscoelasticity.reuss_alv_ortho Function

reuss_alv_ortho(o_phases, fractions) -> NTuple{12, Matrix}

Ortho-form Reuss bound: invert per-phase compliances in ortho form, average, invert back.

MeanFieldHomogenization.Viscoelasticity.dilute_alv_ortho Function

dilute_alv_ortho(o_0, contribs_ortho, fractions) -> NTuple{12, Matrix}

Ortho-form Dilute scheme: o_eff = o_0 + Σ_r f_r · Ñ_r.

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_ortho Function

dilute_dual_alv_ortho(o_0, contribs_compliance_ortho, fractions)
    -> NTuple{12, Matrix}

Ortho-form DiluteDual: invert to compliance ortho form, average, invert back.

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_ortho Function

mori_tanaka_alv_ortho(o_0, A_duts_ortho, contribs_ortho, fractions, f_M)
    -> NTuple{12, Matrix}

Ortho-form Mori-Tanaka: C̃_eff = C̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}, all in ortho form.

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_ortho Function

maxwell_alv_ortho(o_0, contribs_ortho, fractions, o_H_0) -> NTuple{12, Matrix}

Ortho-form Maxwell scheme.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_ortho Function

dilute_concentration_alv_ortho(o_E, o_0, o_P) -> NTuple{12, Matrix}

Ortho-form dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol}.

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_ortho Function

dilute_contribution_alv_ortho(o_E, o_0, o_P) -> NTuple{12, Matrix}

Ortho-form dilute contribution Ñ = ΔC̃ ∘ Ã^dil.

Cracks

MeanFieldHomogenization.Viscoelasticity.cod_kernel_alv Function

cod_kernel_alv(crack::EllipticCrack, C_M_law::ViscoLaw, times;
               Rn = nothing, Rt = nothing) -> NamedTuple

Discrete ALV COD-tensor data for a penny crack η = 1 in an isotropic ALV matrix. Returns the named tuple (B_n = …, B_t = …) of two n × n scalar Volterra matrices (n = length(times)). Each B̃[i,j] approximates the COD coefficient at the time pair (t_i, t_j).

Interface stiffness (Sevostianov-style spring-like interface)

When the crack carries finite interface stiffness kernels Rn(t,t') (normal) and Rt(t,t') (tangential), pass them as scalar ViscoLaws through the Rn / Rt keyword arguments. The traction-free COD matrices B̃_n, B̃_t are then post-corrected via the algebraic identity

B̃_eff = (b · K + B̃^{-1})^{-vol} = B̃ ∘ (𝟙 + b · K ∘ B̃)^{-vol}

(see [69], [59]) where b = semi_minor of the elliptic crack. Limits :

• Rn / Rt = nothing (default) → traction-free penny limit, recovers the existing B_n, B_t.

• Rn, Rt → ∞ (rigid bonding) → B̃_eff_n, B̃_eff_t → 0 (no opening).

Throws if the matrix law is not iso or the crack is not a penny.

MeanFieldHomogenization.Viscoelasticity.compliance_contribution_alv Function

compliance_contribution_alv(crack, C_M_law::ViscoLaw, times) -> Matrix{T}

Discrete (6n × 6n) size-independent compliance contribution H̃ of a penny crack in an isotropic ALV matrix. Computed via the time-space decoupling formula

H̃ = (3/4) · B̃_n · W₁(n̂) + (3/8) · B̃_t · W₆(n̂)

where W₁, W₆ are the canonical Walpole basis tensors of the crack normal axis. When n̂ = e_3 the result is in TI form and routes through the existing TI ALV fast path; arbitrary orientation requires a 6×6 Mandel rotation per (i, j) block (not yet implemented).

Convention: same as the elastic compliance_contribution — the Budiansky-O'Connell density factor is applied separately via delta_compliance_alv.

MeanFieldHomogenization.Viscoelasticity.delta_compliance_alv Function

delta_compliance_alv(crack, H̃, ε) -> Matrix

Apply the Budiansky-O'Connell crack density factor to the size-independent compliance contribution H̃ produced by compliance_contribution_alv, giving the fractional compliance correction ΔJ̃ = (4π/3) ε³ᵈ · H̃ (penny / elliptic geometry, ε³ᵈ = N a b²) — same pre-factor as the elastic case.

MeanFieldHomogenization.Viscoelasticity.stiffness_contribution_alv Function

stiffness_contribution_alv(sphere, C0_law, times) -> Matrix{T}

(6n × 6n) size-independent ALV stiffness contribution of a layered sphere relative to its iso ALV matrix C0_law. Iso parameters (α-, β-blocks of the assembled matrix) are α = 3 Σ_k f_k (M_κ_k − M_κ_0) ∘ α_k, β = 2 Σ_k f_k (M_μ_k − M_μ_0) ∘ β_k, where α_k, β_k are the per-layer localization matrices and M_κ_k, M_μ_k the per-layer Volterra moduli.

The dilute-scheme effective stiffness with this inclusion at volume fraction f is C̃_eff = C̃_0 + f · stiffness_contribution_alv(sphere, …).

stiffness_contribution_alv(crack, C_ref, times) -> Matrix{T}

Discrete (6n × 6n) crack stiffness contribution Ñ = − C̃_ref ∘ H̃ ∘ C̃_ref, mirror of the elastic [stiffness_contribution(crack, C₀)] formula. C_ref may be a ViscoLaw (the matrix law — relaxation auto-built through _trapezoidal_relaxation) or a pre-discretized (6n × 6n) reference matrix (used by SC iterations against the running estimate C_n).

MeanFieldHomogenization.Viscoelasticity.stiffness_contribution_alv_at Function

stiffness_contribution_alv_at(sphere::LayeredSphere, C_ref::AbstractMatrix, times)

Layered-sphere counterpart of stiffness_contribution_alv_at for cracks: the same size-independent ALV contribution, but against a pre-discretized (6n × 6n) reference — the running effective medium of the differential (or self-consistent) ODE rather than the matrix law.

The reference must be isotropic, as everywhere in the layered-sphere ALV recurrences.

stiffness_contribution_alv_at(crack, C_ref::AbstractMatrix) -> Matrix

Variant that takes a pre-discretized (6n × 6n) reference matrix. The compliance contribution is recomputed from the iso parameters of C_ref (only iso ALV matrices are currently supported by compliance_contribution_alv).

MeanFieldHomogenization.Viscoelasticity.delta_stiffness_alv Function

delta_stiffness_alv(crack, Ñ, ε) -> Matrix

Apply the Budiansky-O'Connell crack density factor to the size-independent stiffness contribution Ñ produced by stiffness_contribution_alv, giving the dilute stiffness correction ΔC̃ = (4π/3) ε³ᵈ · Ñ (penny / elliptic) or π ε²ᵈ · Ñ (ribbon). Same pre-factors as the elastic case.

Layered spheres

MeanFieldHomogenization.Viscoelasticity.bulk_localization_alv Function

bulk_localization_alv(sphere, C0_law, times) -> NTuple{N, Matrix}

Per-layer bulk localization matrices α_k(t,t') (n × n each), such that <ε_v>_layer_k = ⟨α_k⟩ ∘ ε_v_∞ in the Volterra sense.

For the bulk Y₀ harmonic, the volume-averaged volumetric strain in layer k is exactly 3 A_k (the mode-2 amplitude B_k has u_r = B_k/r² ⇒ traceless contribution). Therefore α_k = A_k · A_M^{-1} with A_M the matrix-side mode-1 amplitude.

Implementation : amplitude-space recurrence with closed-form interface transitions (bulk_amplitude_seq_alv). Stable for arbitrary modulus contrasts (pores, step-activated layers …).

Reference : Hervé-Zaoui 1993 (elastic) ; ECHOES manual ch07 §"n-layer ALV bulk recurrence" ; sanahuja2013.

MeanFieldHomogenization.Viscoelasticity.bulk_state_seq_alv Function

bulk_state_seq_alv(sphere, C0_law, times)
    -> (inside_states::NTuple{N, Matrix}, s_matrix::Matrix)

Backwards-compatible wrapper that exposes the σ-state form of the bulk recurrence: at each layer k, reconstruct the 2n × n block (u_r ; σ_rr) from the per-layer amplitudes via u_r = A_k r_k + B_k / r_k² and σ_rr = 3 κ_k A_k - 4 μ_k B_k / r_k³.

Most users should rely on bulk_localization_alv directly. This helper is retained for the existing test that asserts the matrix-side state matches the analytical Hervé-Zaoui boundary condition.

MeanFieldHomogenization.Viscoelasticity.bulk_amplitude_seq_alv Function

bulk_amplitude_seq_alv(sphere, C0_law, times)
    -> (NTuple{N, (A_k, B_k)}, A_M, B_M)

Forward-propagate the bulk mode amplitudes (A_k, B_k) (each an n × n Volterra matrix) layer by layer using closed-form interface transitions. Core regularity sets B_1 = 0 and the amplitude seed is A_1 = I_n. Returns the inner-layer amplitudes plus the matrix-side (A_M, B_M).

Supports PerfectInterface, SpringInterface (primal) and MembraneInterface (dual) — same set as the σ-state recurrence — but with no matrix inversion of a fundamental M(r; κ, μ). The only inversions performed are volterra_divide(_, 3κ + 4μ) which remain stable for any non-degenerate modulus.

MeanFieldHomogenization.Viscoelasticity.shear_localization_alv Function

shear_localization_alv(sphere, C0_law, times) -> NTuple{N, Matrix}

Per-layer deviatoric ALV localization matrices β_k(t,t') (n × n each), defined by <ε_d>_layer_k = ⟨β_k⟩ ∘ ε_d_∞ in the Volterra sense — the n × n Volterra matrix that maps a unit deviatoric remote strain to the volume-averaged deviatoric strain in layer k.

The Y₂-harmonic recurrence uses a (4n × 4n) time-major fundamental matrix per layer, two probe states with seed (a, b) = (I, 0) and (0, I) propagated outward, and a final (2n × 2n) Volterra solve that picks the linear combination matching unit far-field (a_{N+1}, b_{N+1}) = (I, 0). Per-layer β_k is the mode-1 amplitude block extracted from the combined inside state.

Reference : ECHOES manual ch07 §"n-layer ALV shear recurrence" ; Hervé-Zaoui 1993 generalized to ALV via [60].

MeanFieldHomogenization.Viscoelasticity.strain_strain_loc_alv Function

strain_strain_loc_alv(sphere, C0_law, times) -> Matrix{T}

(6n × 6n) block matrix describing the volume-averaged strain-strain localization across the entire layered sphere under a unit Volterra far-field strain. In iso form this is A_avg = ⟨α⟩ 𝕁 + ⟨β⟩ 𝕂 with ⟨α⟩ = Σ_k f_k α_k(t,t') and ⟨β⟩ = Σ_k f_k β_k(t,t') (Volterra products).

This is the analog used by the ALV dilute / MT / Maxwell schemes when the inclusion phase is a LayeredSphere.

Order-2 ALV (conductivity / diffusion)

MeanFieldHomogenization.Viscoelasticity.hill_kernel_order2 Function

hill_kernel_order2(ell, K_0_law::ViscoLaw, times) -> Matrix

Build the discrete Hill kernel P̃ (size (3n × 3n)) for an ellipsoidal inclusion in an isotropic ALV matrix. Uses the time-space decoupling formula P̃[block(i,j)] = α₀^{-vol}[i,j] · 𝐈^A.

MeanFieldHomogenization.Viscoelasticity.voigt_alv_order2 Function

voigt_alv_order2(matrices, fractions) -> Matrix

Order-2 Voigt bound: K̃_eff = Σ_r f_r K̃_r.

MeanFieldHomogenization.Viscoelasticity.reuss_alv_order2 Function

reuss_alv_order2(matrices, fractions) -> Matrix

Order-2 Reuss bound: invert each compliance, average, invert back.

MeanFieldHomogenization.Viscoelasticity.dilute_alv_order2 Function

dilute_alv_order2(K_0, contribs, fractions) -> Matrix

Order-2 Dilute scheme: K̃_eff = K̃_0 + Σ_r f_r Ñ_r.

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_order2 Function

dilute_dual_alv_order2(K_0, contribs_compliance, fractions) -> Matrix

Order-2 DiluteDual: invert to compliance space, average, invert back.

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_order2 Function

mori_tanaka_alv_order2(K_0, A_duts, contribs, fractions, f_M) -> Matrix

Order-2 Mori-Tanaka: K̃_eff = K̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}.

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_order2 Function

maxwell_alv_order2(K_0, contribs, fractions; H_0) -> Matrix

Order-2 Maxwell scheme. H_0 is the Hill kernel of the (matrix-only) distribution shape — defaults to a sphere when not specified by homogenize_alv_order2.

MeanFieldHomogenization.Viscoelasticity.homogenize_alv_order2 Function

homogenize_alv_order2(rve, scheme, prop; times)

Backwards-compatible alias for homogenize_alv when the matrix property is order-2. New code should call homogenize_alv directly — the dispatch on order-2 vs order-4 is automatic from the law sample.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_order2 Function

dilute_concentration_alv_order2(K_E, K_0, P) -> Matrix

Order-2 dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔK̃)^{-vol}, all matrices (3n × 3n).

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_order2 Function

dilute_contribution_alv_order2(K_E, K_0, P) -> Matrix

Order-2 dilute contribution Ñ = ΔK̃ ∘ Ã^dil.

MeanFieldHomogenization.Viscoelasticity.iso_order2_params_from_blocks Function

iso_order2_params_from_blocks(M) -> α::Matrix

Decompose a (3n × 3n) block matrix whose every 3×3 block is α[i,j]·𝐈 (iso 2-tensor) into the scalar n × n Volterra matrix α.

MeanFieldHomogenization.Viscoelasticity.iso_order2_blocks_from_params Function

iso_order2_blocks_from_params(α::AbstractMatrix) -> Matrix

Inverse of iso_order2_params_from_blocks: build a (3n × 3n) block matrix α[i,j]·𝐈 per block.

Structured ALV kernel types

MeanFieldHomogenization.Viscoelasticity.AbstractALVKernel Type

AbstractALVKernel{T} <: AbstractMatrix{T}

Abstract supertype for structured ALV kernel wrappers. Concrete subtypes (ALVKernelISO, ALVKernelTI, ALVKernelOrtho) store the symmetry-class parameters compactly and present a (6n × 6n) AbstractMatrix view via lazy getindex.

MeanFieldHomogenization.Viscoelasticity.ALVKernelISO Type

ALVKernelISO{T}(α::Matrix{T}, β::Matrix{T})
ALVKernelISO(M::AbstractMatrix)

Iso ALV kernel: stores the two n × n Volterra parameter matrices α and β (with α = 3K(t,t') and β = 2μ(t,t') in Mandel form). Materializes as a (6n × 6n) block matrix on demand via Matrix(K) or getindex.

MeanFieldHomogenization.Viscoelasticity.ALVKernelTI Type

ALVKernelTI{T}(ℓ::NTuple{6, Matrix{T}}; axis = (0, 0, 1))
ALVKernelTI(M::AbstractMatrix; axis = (0, 0, 1))

TI ALV kernel: stores the six n × n Walpole parameter matrices (ℓ₁, ℓ₂, ℓ₃, ℓ₄, ℓ₅, ℓ₆) with the canonical axis (currently only e₃ supported). Materializes as a (6n × 6n) block matrix on demand.

MeanFieldHomogenization.Viscoelasticity.ALVKernelOrtho Type

ALVKernelOrtho{T}(o::NTuple{12, Matrix{T}}; axes = canonical)
ALVKernelOrtho(M::AbstractMatrix; axes = canonical)

Ortho ALV kernel: stores the twelve n × n parameter matrices of the ortho closure (9 entries of the full unsymmetric 3×3 normal block in Mandel form + 3 shears) with the canonical material frame (e₁, e₂, e₃). Materializes as a (6n × 6n) block matrix on demand.

Orientation symmetrization

MeanFieldHomogenization.Viscoelasticity._maybe_symmetrize_alv Function

_maybe_symmetrize_alv(M, sym) -> Matrix

Apply the orientation-averaging projection corresponding to sym to a (6n × 6n) ALV block matrix : NoSymmetrize (passthrough), IsoSymmetrize (block-by-block exact SO(3) average) and TISymmetrize(axis) (block-by-block exact azimuthal average about the axis).

MeanFieldHomogenization.Viscoelasticity._maybe_symmetrize_alv2 Function

_maybe_symmetrize_alv2(M, sym) -> Matrix

Order-2 counterpart of _maybe_symmetrize_alv, for (3n × 3n) ALV matrices whose blocks are 2-tensors.

Public dispatcher

MeanFieldHomogenization.Viscoelasticity.homogenize_alv Function

homogenize_alv(lam::Laminate, ::Laminated, prop; times, kw...)

Ageing-viscoelastic effective operator of a periodic multilayer, dispatching on the order of the layer laws exactly as the elastic Laminated scheme dispatches on the order of the layer tensors.

homogenize_alv(rve, scheme, prop::Symbol; times) -> Matrix

ALV pipeline: build the discrete block matrices of every phase, compute the ALV Hill kernel and dilute concentration tensors, and dispatch on scheme to the corresponding _alv scheme function.

The function dispatches on the order of the matrix property (read once from the matrix ViscoLaw sample type):

• order-4 (4-tensor / 6×6 Mandel kernel) → returns (6n × 6n) relaxation matrix following the standard Hill-kernel + dilute-concentration pipeline.

• order-2 (2-tensor / 3×3 kernel; conductivity / diffusion / permittivity) → returns (3n × 3n) matrix via the order-2 ALV pipeline (hill_kernel_order2, time-space decoupling).

Supports two inclusion-geometry families (order-4 only):

• Single-shape ellipsoidal (Ellipsoid, Spheroid): the standard Hill-kernel + dilute-concentration pipeline.

• LayeredSphere: bulk + shear ALV recurrences (see bulk_localization_alv and shear_localization_alv) feed stiffness_contribution_alv and strain_strain_loc_alv directly — no Hill kernel is needed. In this case phase_property(rve, name, :C) is ignored (the per-layer moduli stored in the geometry are used instead); pass any ViscoLaw (e.g. heaviside_law(C_0)) as a placeholder.

symmetrize is honored in both orders, with the projector of the matching tensor order (6×6 Mandel blocks for order 4, 2-tensor blocks for order 2 — see _maybe_symmetrize_alv and _maybe_symmetrize_alv2). Voigt and Reuss average the phase matrices directly, so the geometry — and hence symmetrize — does not enter them.

Reference-updating schemes need an isotropic running medium

Every ALV Hill kernel, order 2 and order 4 alike, is built for an isotropic reference: that is what decouples the time and space parts. Voigt, Reuss, Dilute, DiluteDual, MoriTanaka, Maxwell and PCW evaluate it against the fixed, isotropic matrix and accept any inclusion shape or orientation. SelfConsistent and DifferentialScheme evaluate it against their running estimate, so every inclusion phase must keep that estimate isotropic — either a spherical inclusion with an isotropic phase law (LayeredSphere qualifies, its contribution being isotropic by construction), or an isotropic orientation average symmetrize = :iso, which is also what randomly oriented inclusions and cracks mean physically.

An RVE satisfying neither raises an ArgumentError naming the phase, rather than reading iso parameters off a matrix that is not isotropic. A non-isotropic ALV matrix is refused for the same reason. The elastic pipeline has no such restriction — its Hill tensor is available for anisotropic references.

MeanFieldHomogenization.Viscoelasticity.has_visco_property Function

has_visco_property(cell, prop::Symbol = :C) -> Bool

Return true if any member (phase or layer) of cell carries a ViscoLaw under the key prop.

Reads through the raw accessor cell_container_property: this is a type inspection, and going through the resolving accessor would run a whole inner homogenization on every declaratively nested cell just to look at a type. A nested cell is instead recursed into, so a ViscoLaw buried one scale down is still reported.