API — LayeredSpheroid
layer_count, layer_modulus, layer_interface and layer_volume_fraction are shared generics extended from LayeredSpheres — see API — LayeredSphere for their docstrings; they apply unchanged to LayeredSpheroid.
So are the pointwise-field generics get_layer, local_temperature, local_gradient, local_flux and the four local_*_*_loc couplings: one binding carries both inclusion families, and all their docstrings — the sphere's methods and the spheroid's — are rendered together on that same page.
MeanFieldHomogenization.LayeredSpheroids Module
MeanFieldHomogenization.LayeredSpheroidsIsotropic n-layer confocal spheroidal composite inclusion (core + concentric confocal shells), in conduction (thermal / electric / Darcy) and in axisymmetric elasticity (spheroid_elastic_coefficients, case I of Duan et al. 2005, with perfect interfaces).
The two do not share a solver shape. Confocal surfaces are not homothetic, so the harmonic degrees couple — in conduction only at an imperfect interface, which is why its blocks are diagonal in the perfect case, but in elasticity at every interface. There is no per-interface transfer matrix to chain on the elastic side; elasticity.jl assembles one global system instead. Public entry points: LayeredSpheroid, layered_spheroid_from_fractions.
Like LayeredSphere, a composite spheroid has no Hill tensor — it plugs into the mean-field schemes through its volume-averaged concentration (gradient/flux) tensors, assembled layer by layer via a confocal spheroidal-harmonic transfer-matrix recurrence (conductivity.jl) instead of the sphere's simple 2×2 state-vector propagation. Perfect (PerfectInterface), Kapitza (KapitzaInterface, resistance) and surface-conductive (SurfaceConductiveInterface, conductance) interfaces — reused from LayeredSpheres — couple different harmonic degrees, unlike the sphere, requiring the truncated series machinery of legendre.jl / coupling.jl.
That path into the schemes now works in elasticity too, prolate and oblate alike: spheroid_strain_concentration assembles the full transversely isotropic concentration tensor from the three elementary problems, and strain_strain_loc / stiffness_contribution hand it to the schemes. Only perfect interfaces, and only the default axis ê₃.
MeanFieldHomogenization.LayeredSpheroids.LayeredSpheroid Type
LayeredSpheroid{T, N, Q, Cs, Is} <: AbstractLayeredInclusion{3, T}Isotropic N-layer confocal spheroidal composite inclusion (core + concentric confocal shells). Carries the geometry for both physics — conduction and axisymmetric elasticity; the moduli tuple holds TensISO{2,3} conductivities or TensISO{4,3} stiffnesses accordingly. Type parameters:
T— element type of the geometric scalars (radii, focal distance).N— number of layers (≥ 1).Q— element type of the confocal parameterq:Tfor prolate,Complex{T}for oblate (the complex substitutionc → -i c̄,q → i τ).Cs,Is— concrete types of themoduli/interfacestuples.
Use the keyword constructor LayeredSpheroid(axis_radii, disk_radii, moduli; interfaces, Nseries, axis) for the common case, or layered_spheroid_from_fractions to specify layers by volume fraction.
MeanFieldHomogenization.LayeredSpheroids.layered_spheroid_from_fractions Function
layered_spheroid_from_fractions(ω, outer_axis_radius, layer_fractions, moduli;
interfaces, Nseries = 5, axis = (0., 0., 1.))Convenience constructor specifying layers by volume fraction: build an N-layer confocal spheroid of given outer aspect ratio ω (> 1 prolate, < 1 oblate) and outer axis semi-axis outer_axis_radius, with each layer occupying the prescribed fraction of the total volume (layer_fractions, normalized to sum 1, core first).
The geometry itself comes from confocal_layer_radii, so the same semi-axes are available to build a finite-element counterpart of this inclusion.
MeanFieldHomogenization.LayeredSpheroids.confocal_layer_radii Function
confocal_layer_radii(ω, outer_axis_radius, layer_fractions) -> (axis_radii, disk_radii)The per-layer semi-axes of an N-layer confocal spheroid of outer aspect ratio ω (> 1 prolate, < 1 oblate) and outer axis semi-axis outer_axis_radius, each layer occupying the prescribed fraction of the total volume (core first, normalized to sum 1). Inner confocal parameters come from bisection on the volume relation V(q) ∝ |q(q²-1)|.
Returned in the argument order LayeredSpheroid takes — axis (revolution) semi-axes first, transverse second — so the same pair builds the analytic inclusion and a finite-element one, and the two cannot describe different bodies. That is the whole reason this is a function rather than a step inside a constructor: a finite-element cross-check of the analytic solution is only worth anything if both are handed the same geometry.
MeanFieldHomogenization.LayeredSpheroids.layer_q Function
layer_q(spheroid, k) -> QConfocal parameter q_k of the outer boundary of layer k.
MeanFieldHomogenization.LayeredSpheroids.layer_semiaxes Function
layer_semiaxes(spheroid, k) -> (axis, disk)Real (axis, disk) semi-axes of the outer boundary of layer k.
MeanFieldHomogenization.LayeredSpheroids.outer_semiaxes Function
outer_semiaxes(spheroid) -> (axis, disk)Real (axis, disk) semi-axes of the outermost boundary q_N.
MeanFieldHomogenization.LayeredSpheroids.spheroid_state_sequence Function
spheroid_state_sequence(s, k₀, trans) -> Vector{Vector}Series coefficient vectors Xℓ = [Aℓ; Bℓ] (eq:AlBlXl), ℓ = 1, …, N+1 (the (N+1)-th being the matrix), for the axial (trans = false) or transverse (trans = true) problem in an isotropic matrix k₀. B₁ = 0 (core regularity, eq:axiBCbi) and A_{N+1} = (±1, 0, …, 0) (unit remote field, + axial / − transverse, eq:axiBCai/eq:transBCai) are imposed exactly.
MeanFieldHomogenization.LayeredSpheroids.spheroid_ba_ratios Function
spheroid_ba_ratios(s, k₀) -> (ba_axial, ba_trans)The two ratios b^0_{N+1,1}/a^0_{N+1,1} and b^1_{N+1,1}/a^1_{N+1,1} (eq:axiasb/eq:transasb) driving the volume-averaged concentration tensors. Kept in their native (possibly complex, for the oblate substitution q = iτ) type — casting to real is only valid on the FINAL shape-function product (ba · 𝒯/𝒰(q_N), done in spheroid_gradient_gradient / spheroid_flux_gradient / _spheroid_concentration), never on ba alone.
MeanFieldHomogenization.LayeredSpheroids.LayeredSpheroidTransportFields Type
LayeredSpheroidTransportFields(s, k₀)Precomputed pointwise transport solution of a LayeredSpheroid in the isotropic matrix k₀: the axial and transverse confocal-harmonic coefficient sequences, solved once.
Build it once and pass it wherever (s, k₀) is accepted. The (s, k₀) forms re-run both state sequences on every call, which a field map over thousands of points pays for in full.
Transport twin of LayeredSphereTransportFields.
MeanFieldHomogenization.LayeredSpheroids.coupling_matrices Function
coupling_matrices(q, Nseries; method = :quadrature) -> (I, J, K, L)The four Nseries × Nseries interface coupling matrices at confocal parameter q (real for prolate, q = iτ for oblate), restricted to the odd degrees 1, 3, …, 2·Nseries-1 (index r ↔ degree 2r-1).
method = :quadrature (default) integrates the paper's closed-form definitions directly (stable in Float64); method = :series uses the BigFloat monomial-coefficient summation of the original implementation (see the module docstring).
MeanFieldHomogenization.LayeredSpheroids.legendre_odd Function
legendre_odd(kind::Symbol, x, Nseries::Int) -> (vals, derivs)Values and derivatives of the requested Legendre kind at the Nseries ODD degrees 1, 3, …, 2·Nseries − 1, as length-Nseries Vectors (index r ↔ degree 2r − 1).
kind ∈ (:P0, :Q0, :P1, :P1p, :Q1):
:P0—Pₙ(x)(m=0, any branch):Q0—Qₙ(x)(m=0, q branch, |x|>1):P1—Pₙ¹(x)(m=1, q branch, |x|>1):P1p—Pₙ¹(x)(m=1, p branch, |x|≤1):Q1—Qₙ¹(x)(m=1, q branch, |x|>1)
MeanFieldHomogenization.LayeredSpheroids.legendre_table Function
legendre_table(kind::Symbol, x, Nmax::Int) -> (vals, derivs)Values and derivatives of the requested Legendre kind at every degree 0, 1, …, Nmax, as length-Nmax + 1 Vectors (index n + 1 ↔ degree n).
kind ∈ (:P0, :Q0, :P1, :P1p, :Q1):
:P0—Pₙ(x)(m = 0, any branch):Q0—Qₙ(x)(m = 0, q branch,|x| > 1):P1—Pₙ¹(x)(m = 1, q branch,|x| > 1):P1p—Pₙ¹(x)(m = 1, p branch,|x| ≤ 1):Q1—Qₙ¹(x)(m = 1, q branch,|x| > 1):P2—Pₙ²(x)(m = 2, q branch,|x| > 1):P2p—Pₙ²(x)(m = 2, p branch,|x| ≤ 1):Q2—Qₙ²(x)(m = 2, q branch,|x| > 1)
MeanFieldHomogenization.LayeredSpheroids.legendre_degrees Function
legendre_degrees(kind::Symbol, x, degrees) -> (vals, derivs)Values and derivatives at an arbitrary set of degrees, in the order given.
Conduction needs the odd degrees alone, which is what legendre_odd selects. Elasticity does not: the case-I Papkovich–Neuber potentials split by parity — φ₀ on the even degrees, φ₃ on the odd ones — so the degree list is part of the problem rather than a property of the module. The underlying tables are built for 0:maximum(degrees) either way, the recurrences being what they are, so asking for a sparse set costs nothing extra.
Qₙ depends on the highest degree requested
Pₙ grows upward and is bit-identical whatever Nmax is asked for. Qₙ is the minimal solution and is built by Miller's downward recurrence, which starts above the highest degree requested and normalizes on a closed form — so the arithmetic path, and with it the last couple of bits, depends on Nmax. Two calls with different degree sets agree to a few ulp, never bit-for-bit.
MeanFieldHomogenization.LayeredSpheroids.AxisymmetricCase Type
AxisymmetricCase()Case I of Duan et al. (2005): a remote strain diag(ε_t, ε_t, ε_a) about the spheroid's axis. The Papkovich–Neuber gauge is φ₁ = φ₂ = 0, both surviving potentials are of order m = 0, and the fields are axisymmetric — u_φ and σ_φq vanish structurally rather than by assumption.
A singleton, so the two remaining elementary problems (transverse and longitudinal shear, which need orders m = 2 and m = 1) can be added by dispatch without disturbing this one.
MeanFieldHomogenization.LayeredSpheroids.spheroid_elastic_coefficients Function
spheroid_elastic_coefficients(s, C₀, εa, εt; case = AxisymmetricCase(),
ngauss = 0) -> NamedTupleSolve the elastic n-layer confocal spheroid under an axisymmetric remote strain diag(εt, εt, εa) about the spheroid's axis.
Returns (; modes, amplitudes, residual): modes[ℓ] lists the (kind, degree) pairs of region ℓ (1:N the layers, N+1 the matrix) and amplitudes[ℓ] their values, with the matrix's prescribed remote part left out. residual is the relative least-squares residual, which is a diagnostic and not a fitting error — see the note below.
Why the system is solved in least squares
Each interface contributes 4𝒩 conditions and each region 4𝒩-1 amplitudes, degree 0 of φ₀ being a constant potential that moves nothing. The assembled system is therefore over-determined by one row per interface, and those rows are redundant, not conflicting: the residual comes out at machine precision. Solving in least squares is the honest way to use that, and a residual that stops being negligible is a signal that something upstream is wrong.
How far Nseries is worth pushing
A single homogeneous spheroid is exact at any truncation — its series collapses to two coefficients. A genuinely layered one converges geometrically, the residual falling by roughly a factor 3 per unit of Nseries. In Float64 that stops paying at about Nseries = 12 (relative accuracy ~1e-7): past there the least-squares conditioning dominates, the residual stops falling, and the answer parts company with a 256-bit one. Raise the element type rather than the truncation. The conduction solver records the same limit after Barthélémy & Bignonnet's appendix C; the elastic blocks are wider, so it arrives sooner.
Only PerfectInterface is supported so far; imperfect interfaces need their jump terms added to the four conditions. Oblate spheroids are refused: their confocal parameter is complex and nothing here has been checked against a reference for it.
MeanFieldHomogenization.LayeredSpheroids.spheroid_core_strain Function
spheroid_core_strain(s, C₀, εa, εt; kw...) -> (εa_in, εt_in)Uniform strain in the core of an elastic confocal spheroid under the axisymmetric remote strain diag(εt, εt, εa).
The core carries only the regular modes, and the two that survive — degree 2 of φ₀ and degree 1 of φ₃ — are exactly a uniform strain, read back through the same relations that build the remote field:
ε_t = -3A₂ / (4 c² μ₁), ε_a = -γ(1 - 2ν₁)/μ₁ - 2ε_t, γ = C₁/c.For a single layer this is Eshelby's result, which is what the tests check it against.
MeanFieldHomogenization.LayeredSpheroids.TransverseShearCase Type
TransverseShearCase(trig)Case II: a remote transverse shear. trig = :sin loads 2\varepsilon_{12}, trig = :cos loads \varepsilon_{11}-\varepsilon_{22}; the two are the same problem rotated by
Active potentials are φ₀, φ₃ at order m = 2 and the tied pair (φ₁, φ₂) at m = 1.
MeanFieldHomogenization.LayeredSpheroids.LongitudinalShearCase Type
LongitudinalShearCase(trig)Case III: a remote shear between the axis and a transverse direction. trig = :cos loads 2\varepsilon_{13}, trig = :sin loads 2\varepsilon_{23}.
Active potentials are φ₀, φ₃ at order m = 1 and φ₁ (resp. φ₂) at m = 0. No rigid-body rotation is added — see this file's header.
MeanFieldHomogenization.LayeredSpheroids.ModeGroup Type
ModeGroup(parts, degrees)Modes that share a single amplitude, and the degrees they run over.
parts is a list of (potential, m, trig, weight). Most groups have one part. Case II needs more: its remote field ties φ₁ and φ₂ together — Duan's "
MeanFieldHomogenization.LayeredSpheroids.PNMode Type
PNMode(potential, regularity, n, m, trig)One harmonic mode of one Papkovich–Neuber potential.
potential ∈ 0:3— which ofφ₀, φ₁, φ₂, φ₃it belongs to.regularity ∈ (:regular, :irregular)— theq-dependence, in Barthélémy & Bignonnet's own words (2020, §2):Pₙᵐ(p)Pₙᵐ(q)has a finite limit asq → 1and is a regular harmonic,Pₙᵐ(p)Qₙᵐ(q)blows up there and is an irregular one. A core admits only regular harmonics, the matrix only irregular ones plus the remote field.n,m— degree and order.trig ∈ (:cos, :sin)— the azimuthal factorcos mφorsin mφ.
Amplitude naming
Barthélémy & Bignonnet write the series of layer ℓ as
Σ Pₙᵐ(p) { [aᵐ_{ℓ,n} Pₙᵐ(q) + bᵐ_{ℓ,n} Qₙᵐ(q)] cos mφ
+ [cᵐ_{ℓ,n} Pₙᵐ(q) + dᵐ_{ℓ,n} Qₙᵐ(q)] sin mφ }so a, b, c, d mean regular-cos, irregular-cos, regular-sin, irregular-sin — see bb_letter. Conduction has one field; elasticity has four potentials, so the letters carry a potential index as well, a^{i,m}_{ℓ,n} … d^{i,m}_{ℓ,n}.
A letter encodes (regularity, trig), never a potential
An earlier version of this module labeled case I's four families :A, :B, :C, :D for (φ₀, reg), (φ₀, irr), (φ₃, reg), (φ₃, irr). That collides with the convention above, where c and d are the sine families. The pair is carried explicitly here for that reason.
MeanFieldHomogenization.LayeredSpheroids.bb_letter Function
bb_letter(mode) -> SymbolBarthélémy & Bignonnet's amplitude letter for a mode: :a regular-cos, :b irregular-cos, :c regular-sin, :d irregular-sin. The potential index rides alongside as a superscript and is not part of the letter.
MeanFieldHomogenization.LayeredSpheroids.mode_fields Function
mode_fields(mode, ϕ, p, q, c, μ, ν, ::Type{T}) -> (u, t)Displacement and traction of one unit-amplitude mode, in the chart frame: u = (u_φ, u_p, u_q) and t = σ·e_q = (σ_φq, σ_pq, σ_qq).
This is the single evaluator all three elementary problems go through. What distinguishes them is which modes are in the list, not how a mode is turned into fields.
MeanFieldHomogenization.LayeredSpheroids.spheroid_strain_concentration Function
spheroid_strain_concentration(s::LayeredSpheroid, C₀; D = 6, kw...) -> Tens{4,3}Volume-averaged strain concentration tensor 𝔸 of an elastic n-layer confocal spheroid in an infinite isotropic matrix C₀, defined by ⟨ε⟩ = 𝔸 : E over the whole composite inclusion.
This is what a mean-field scheme consumes. It is assembled from six solves, one per independent remote strain — _basis_loadings lists them — each averaged by the surface integral of _avg_strain. The result is returned as a general fourth-order tensor rather than a TensTI{4}: it is transversely isotropic about the spheroid's axis, but the tests check that rather than the type asserting it, and 𝔸 has no major symmetry to exploit.
Returns (; A, residuals), residuals being the six least-squares residuals — diagnostics, machine-precision for a single inclusion.
The axis
𝔸 comes out in the global frame, so a spheroid whose axis is not ê₃ is not yet handled here: the elementary problems are written about ê₃. Build the spheroid with the default axis and rotate the result.
MeanFieldHomogenization.LayeredSpheroids.spheroid_layer_strain_concentration Function
spheroid_layer_strain_concentration(s, C₀; D = 6, kw...) -> (; A, layers, residuals)Per-layer strain concentration tensors of an elastic confocal spheroid, plus their volume-weighted sum.
layers[k] is 𝔸_k with ⟨ε⟩_k = 𝔸_k : E averaged over layer k alone, and A is Σ_k f_k 𝔸_k. The two are computed independently — the layers from differences of surface moments on their own boundaries, A from a single moment on the outer one — so their agreement is a check rather than a restatement.
Per-layer averages are what a stiffness contribution tensor needs: N_C = Σ_k f_k (ℂ_k - ℂ_0) : 𝔸_k cannot be recovered from the total alone, the layers having different moduli.