API — Core

MeanFieldHomogenization Module

MeanFieldHomogenization

Julia package for mean-field homogenization of heterogeneous materials.

MeanFieldHomogenization unifies the computation of Hill polarization tensors for ellipsoidal inhomogeneities, crack opening displacement (COD) tensors, stress and displacement intensity factors, homogenization schemes over representative volume elements (RVEs), and ageing viscoelastic constitutive laws, sharing a common abstraction for inclusions, algorithms, and material symmetry classes.

Sub-modules

• MeanFieldHomogenization.Elliptic — type-generic Legendre and Carlson elliptic integrals (ForwardDiff- and Sym-compatible).

• MeanFieldHomogenization.Core — abstractions (AbstractInclusion, AbstractAlgorithm, MaterialSymmetry), shared numerics (Green / Newton kernels, Masson-style residue, DECUHR integrand), modulus extractors, and central dispatch.

• MeanFieldHomogenization.Elasticity — Hill polarization for ellipsoidal inclusions and infinite cylinders (2D / 3D, isotropic and anisotropic matrix).

• MeanFieldHomogenization.Cracks — COD tensors, compliance contributions, SIF and DIF for elliptic and ribbon cracks.

• MeanFieldHomogenization.Conductivity — 2nd-order Hill tensor for conductivity / diffusion problems.

• MeanFieldHomogenization.Interactions — two-inclusion interaction tensors Γ^{ab}, the shared ingredient of the N-body schemes (equivalent inclusion method, cluster model): closed forms for ball and disk pairs, multipole expansion for general ellipsoids, and periodic lattice sums.

• MeanFieldHomogenization.LayeredSpheres — n-layer composite spheres with five interface types, volume-average and pointwise localization.

• MeanFieldHomogenization.LayeredSpheroids — n-layer confocal spheroids in conduction, with imperfect interfaces.

• MeanFieldHomogenization.Schemes — RVEs, amounts, symmetrization and the homogenization schemes themselves (dilute, Mori–Tanaka, self-consistent, PCW, Maxwell, differential).

• MeanFieldHomogenization.Poromechanics — poroelastic upscaling of a saturated medium with a homogeneous solid phase: Biot tensor and modulus, drained ↔ undrained conversion, Skempton tensor, effective stresses. A post-processor of a homogenized stiffness, not a scheme.

• MeanFieldHomogenization.Viscoelasticity — linear viscoelasticity by two routes: the ageing time-domain one through Volterra operators (homogenize_alv), and the non-ageing Laplace-Carson one (homogenize_lc), with a catalog of rheological models, numerical Laplace inversion, and the exact Kelvin ↔ Maxwell conversion joining them.

• MeanFieldHomogenization.CustomInclusions — the user-defined inclusion contract: CustomInclusion and check_inclusion_interface.

• MeanFieldHomogenization.FiniteElements — inclusions whose response comes out of a finite-element resolution of the Eshelby problem (FEEllipticCrack, FEExcenteredSphere); the discretization comes from a backend extension, MeanFieldHomogenizationFerriteExt or MeanFieldHomogenizationGridapExt.

• MeanFieldHomogenization.NeuralInclusions — inclusions whose response comes out of a trained neural network (NeuralHillInclusion, NeuralLocalizationInclusion), together with the sampling and fitting machinery; the optimizer comes from MeanFieldHomogenizationLuxExt, evaluation needs no extra dependency.

• MeanFieldHomogenization.Constitutive — the package as a constitutive law at each Gauss point of a structural finite-element computation, the role an MFront behavior or an Abaqus UMAT plays: material_response, per-point internal state, consistent tangent, and the Tensors.jl bridge. The mirror image of FiniteElements — there the finite elements are inside MFH, here MFH is inside the finite-element code.

Shared generic interface

hill_tensor(ell::AbstractEllipsoidalInclusion, C₀; method=:auto, ...)
cod_tensor(crack::AbstractCrack, C₀; method=:auto, ...)
compliance_contribution(crack, C₀; method=:auto, ...)    # returns H (or R)
delta_compliance(crack, H, ε)                             # ΔS = factor · ε · H
delta_resistivity(crack, R, ε)                            # ΔR = factor · ε · R
sif(crack, C₀, Σ; method=:auto, ...)
dif(crack, C₀, Σ; method=:auto, ...)

All high-level entry points share the same algorithmic traits (Analytical, Residue, DECUHR) and the same material-symmetry dispatch rules. See the developer documentation (docs/src/developer/) for guidance on extending the package with new inclusions, algorithms or schemes.

MeanFieldHomogenization.Core Module

MeanFieldHomogenization.Core

Shared abstractions and numerical kernels used throughout MeanFieldHomogenization.

Contents

• counters.jl : opt-in work counters for the benchmark harness

• abstractions.jl : inclusion hierarchy (AbstractInclusion …)

• cells.jl : the homogenization cell contract (AbstractHomogenizationCell, AbstractParameter) and the declarative multiscale seam (Homogenized, NestedParameter)

• traits.jl : algorithm and material-symmetry traits

• bases.jl : helpers around TensND bases

• tensor_helpers.jl : low-level utilities (_δ, _C_array, Voigt)

• moduli.jl : modulus extractors for the common symmetry classes

• newton_potential.jl : Newton potentials (2D / 3D)

• green_kernel.jl : closed-form 3×3 inverse (_inv3)

• laminate_algebra.jl : in-plane / out-of-plane Kelvin-Mandel block algebra of a periodic laminate (plane_pinv, flat_hill, laminate_stiffness, …)

• green_residue.jl : Masson / Cauchy residue summation

• green_helpers.jl : quadrature-agnostic Green-function helpers

• green_dipole.jl : real-space Kelvin Green gradient and the dipole far field of a polarized inclusion (isotropic)

• green_operator.jl : real-space Green operator 𝔾⁰(x) — one derivative beyond green_dipole.jl, the regular kernel of the Lippmann-Schwinger equation (isotropic)

• green_aniso.jl : its anisotropic counterpart — Barnett line integral for elasticity, closed form for conduction — and the green_operator dispatcher over both

• quadrature.jl : DECUHR cubature backend seam

• dispatch.jl : central _resolve_algo mechanism

• dispatch_pair.jl : _resolve_pair_algo — its two-inclusion counterpart

MeanFieldHomogenization.Core.AbstractInclusion Type

AbstractInclusion{T<:Number}

Root abstract supertype for every inclusion geometry recognized by MeanFieldHomogenization. The type parameter T is the element type of the geometric scalars (semi-axes, half-widths, …) and propagates through every tensor produced by the package, supporting Float64, ForwardDiff.Dual, SymPy.Sym, Symbolics.Num, …

MeanFieldHomogenization.Core.AbstractEllipsoidalInclusion Type

AbstractEllipsoidalInclusion{dim,T} <: AbstractInclusion{T}

Supertype for ellipsoidal inclusions — solid ellipsoids (and their degenerate limits: spheres, cylinders, discs …). The first type parameter dim encodes the spatial dimension (2 or 3).

MeanFieldHomogenization.Core.AbstractCrack Type

AbstractCrack{T} <: AbstractInclusion{T}

Supertype for flat-crack geometries (elliptic / ribbon / penny). Cracks always live in 3D physical space — the spatial dimension of the geometry (= 2 in the crack plane) is not exposed in the type, as downstream algorithms uniformly operate on the 3D stiffness tensor.

MeanFieldHomogenization.Core.AbstractLayeredInclusion Type

AbstractLayeredInclusion{dim,T} <: AbstractInclusion{T}

Supertype for multi-layer inclusions. The concrete LayeredSphere (concentric isotropic shells, Hervé-Zaoui recurrences for bulk, shear and conductivity, with five interface types) is shipped in the LayeredSpheres sub-module and extended to the ageing-viscoelastic setting in Viscoelasticity/layered_alv.jl. Open extensions (coated cylinders, anisotropic layers, excentered spheres) are tracked in docs/src/developer/roadmap.md.

MeanFieldHomogenization.Core.AbstractAlgorithm Type

AbstractAlgorithm

Abstract supertype for computation algorithms exposed through MeanFieldHomogenization. The subtypes are plain singleton structs used as dispatch tags by the internal _kernel(...) machinery.

MeanFieldHomogenization.Core.Auto Type

Placeholder singleton used by _resolve_algo(Val(:auto), …) as an explicit automatic selection request.

MeanFieldHomogenization.Core.Analytical Type

Closed-form / analytical algorithm — used for every isotropic matrix and for every 2D / conductivity case.

MeanFieldHomogenization.Core.Residue Type

Residue-theorem algorithm — polynomial-root based, Float64 only (incompatible with ForwardDiff and SymPy).

MeanFieldHomogenization.Core.DECUHR Type

DECUHR 2D hyper-cubature (Espelid & Genz 1994) via Integrals.jl DecuhrAlgorithm. Suitable for general anisotropy in 3D.

MeanFieldHomogenization.Core.NestedQuadGK Type

Nested 1D QuadGK cubature (fallback to DECUHR, ForwardDiff-compatible). Historically shipped as DECUHR before the split.

MeanFieldHomogenization.Core.CylinderQuadrature Type

1D QuadGK quadrature dedicated to infinite cylinders (transverse-plane parametrization ζ(φ) = (0, cos φ / b, sin φ / c)). ForwardDiff-compatible, selected whenever an AbstractEllipsoidalInclusion with a cylindrical-shape trait meets a general-anisotropic 3D stiffness.

MeanFieldHomogenization.Core.Multipole Type

Multipole <: AbstractAlgorithm

Truncated multipole expansion of a two-inclusion interaction integral (brisard2014, §4.2): the regular part of the Green operator is Taylor-expanded about the line of centers and the resulting monomials are integrated over the two regions in closed form. Exact for balls; an asymptotic series in (size / separation) for a general ellipsoid.

MeanFieldHomogenization.Core.MaterialSymmetry Type

MaterialSymmetry

Abstract supertype classifying the symmetry class of a TensND tensor. Used by dispatch rules to select the most specific algorithm (e.g. closed-form formulas for isotropic and transversely isotropic matrices).

MeanFieldHomogenization.Core.material_symmetry Function

material_symmetry(C₀)

Return the MaterialSymmetry trait corresponding to the TensND tensor C₀. Dispatches on the concrete TensND type: TensISO → IsotropicSym, TensTI{4} → TransverselyIsotropicSym, TensOrtho → OrthotropicSym, anything else → GeneralAnisotropicSym.

MeanFieldHomogenization.Core.IsotropicSym Type

Isotropic material (TensND TensISO).

MeanFieldHomogenization.Core.TransverselyIsotropicSym Type

Transversely isotropic material (TensND TensTI{4} / TensTI).

MeanFieldHomogenization.Core.OrthotropicSym Type

Orthotropic material (TensND TensOrtho).

MeanFieldHomogenization.Core.GeneralAnisotropicSym Type

General anisotropic material (generic Tens subtype with no structured symmetry).

MeanFieldHomogenization.Core.extract_iso_conductivity Function

extract_iso_conductivity(K₀::TensISO{2,d}) -> k

Extract the (scalar) conductivity coefficient of an isotropic 2nd-order transport tensor TensISO{2,d} (  ).

MeanFieldHomogenization.Core.extract_iso_moduli Function

extract_iso_moduli(C₀::TensISO{4,3}) -> (E, ν)

Extract Young's modulus E and Poisson's ratio ν from an isotropic 4th-order stiffness TensISO{4,3}. The internal TensND convention is      i.e. C₀.data = (3k, 2μ).

extract_iso_moduli(C₀::TensISO{4,2}) -> (E, ν)

2D plane-strain counterpart of extract_iso_moduli for TensISO{4,2}. The same 3D formulas are used because the TensND storage is dimension-agnostic (TensISO{4,d} stores the same (α, β) = (3k, 2μ) pair).

MeanFieldHomogenization.Core.extract_ti_moduli Function

extract_ti_moduli(C₀, n̂) -> (E, H, ν₁, ν₂, Γ)

Read the five TI compliance moduli out of a stiffness tensor C₀ whose axis of symmetry is n̂. Used by the closed-form COD formulas of the Cracks sub-module.

The moduli are defined through the compliance tensor   as:

•  

•   

•    

•    

•     

See the package documentation for the full derivation.

MeanFieldHomogenization.Core.newton_potential_3d Function

newton_potential_3d(a, b, c) -> (Iv, IIv)

Newton potential integrals for a 3-D ellipsoid with semi-axes a ≥ b ≥ c > 0.

Returns:

• Iv = (I_a, I_b, I_c) with I_a + I_b + I_c = 4π

• IIv = (I_aa, I_bb, I_cc, I_bc, I_ca, I_ab)

Two methods are provided:

• T<:Real (includes Float64, ForwardDiff.Dual): numerically stable case-split using tolerance comparisons (sphere, oblate, prolate, triaxial with elliptic integrals).

• T<:Number (e.g. SymPy.Sym, Symbolics.Num): structural equality via isequal selects the same four cases.

newton_potential_3d(ell::Ellipsoid{3})

Ellipsoid-level convenience wrapper that forwards to Core.newton_potential_3d.

MeanFieldHomogenization.Core.newton_potential_2d Function

newton_potential_2d(a, b) -> (Ia, Ib)

Newton potential integrals for a 2-D ellipse with semi-axes a ≥ b > 0. Returns (Ia, Ib) with Ia + Ib = 2π.

Works for any T<:Number including ForwardDiff.Dual, SymPy.Sym, Symbolics.Num.

Formulas:

• Circle: Ia = Ib = π

• General ellipse: Ia = 2πb/(a+b), Ib = 2πa/(a+b)

newton_potential_2d(ell::Ellipsoid{2})

Ellipsoid-level convenience wrapper that forwards to Core.newton_potential_2d.

MeanFieldHomogenization.Core.newton_potential_3d_cylinder Function

newton_potential_3d_cylinder(b, c) -> (Iv, IIv)

Newton potential integrals for an infinite cylinder of elliptic cross-section with transverse semi-axes b ≥ c > 0 (cylinder axis = e₁, transverse plane = (e₂, e₃)).

Obtained as the limit a → ∞ of newton_potential_3d — the cylinder axis contributes no finite Newton mass (I_a = I_aa = I_ab = I_ac = 0) and the transverse potentials collapse to simple rational expressions in (b, c).

Returns:

• Iv = (I_a, I_b, I_c) with I_a + I_b + I_c = 4π, I_a = 0.

• IIv = (I_aa, I_bb, I_cc, I_bc, I_ca, I_ab) with I_aa = I_ab = I_ac = 0.

Two methods are provided:

• T<:Real (includes Float64, ForwardDiff.Dual): numerically stable case-split via tolerance comparison to pick the circular (b = c) or the elliptic (b > c) branch.

• T<:Number (SymPy.Sym, Symbolics.Num, …): structural equality via isequal selects the two branches.

Both branches are written as closed-form limits — no 1/(b² − c²) style denominators, so the routine is free of 0/0 indeterminacies at b = c and differentiable through ForwardDiff.

newton_potential_3d_cylinder(cyl::Cylinder)

Cylinder-level convenience wrapper that forwards to Core.newton_potential_3d_cylinder using the transverse semi-axes of cyl.

MeanFieldHomogenization.Core.dimension Function

dimension(incl::AbstractInclusion) -> Int

Spatial dimension of the inclusion's ambient space (2 or 3 for the concrete inclusions shipped with the package).

MeanFieldHomogenization.Core.element_type Function

element_type(incl::AbstractInclusion{T}) -> Type{T}

Element type of the geometric scalars stored in the inclusion (Float64, ForwardDiff.Dual, SymPy.Sym, …).

MeanFieldHomogenization.Core.inclusion_basis Function

inclusion_basis(incl::AbstractInclusion) -> TensND.AbstractBasis

Local principal basis of the inclusion (principal frame for an ellipsoid, for a crack, …). Used by downstream algorithms to rotate the matrix stiffness / conductivity into the inclusion frame.

MeanFieldHomogenization.Core.shape_trait Function

shape_trait(incl::AbstractInclusion) -> Type

Concrete shape classification of the inclusion, used as a type parameter for Holy-style dispatch in the downstream kernels. Typical values: Spherical, Prolate, Oblate, Triaxial, Circular, Elliptic (ellipsoids), Penny, EllipticShape, Ribbon (cracks).

MeanFieldHomogenization.Core.shape_tensor Function

shape_tensor(incl::AbstractInclusion) -> AbstractTens{2}

Symmetric 2nd-order tensor encoding both the semi-axes and the orientation of an equivalent ellipsoidal envelope of the inclusion, in the global (canonical) frame.

Optional

No kernel in the package reads it: an inclusion that supplies its own response tensors owes nothing about its outer shape, and a morphology with no ellipsoidal envelope simply has no shape_tensor. Implement it when the notion is meaningful.

where is the rotation matrix mapping the canonical frame onto the inclusion's local basis and the diagonal entries are the semi-axes in the order dictated by the local basis.

Conventions for degenerate cases:

Inclusion	Diagonal (principal frame)
Ellipsoid{3}	(a₁, a₂, a₃)
Ellipsoid{2}	(a₁, a₂)
Cylinder	(Inf, b, c) — axis
EllipticCrack	(a, b, 0) — normal
RibbonCrack	(Inf, b, 0)

Exact rotation-group averages

Exact averages of a tensor over a rotation group — the exact counterpart of the best-fit projections in API — Schemes. See Symmetrization showcase for the comparison between the two.

These live in TensND — they are pure tensor algebra, with nothing homogenization-specific about them — and are re-exported by this package, so MeanFieldHomogenization.isotropify and Core.isotropify both resolve to them:

Name	What it returns
isotropify(t)	the exact SO(3) average, TensISO{4} or TensISO{2}
transverse_isotropify(t, n)	the exact azimuthal average about n, TensTI{4,T,8} or TensTI{2,T,3} — ℓ₃ ≠ ℓ₄ and the antisymmetric couplings ℓ₇, ℓ₈ are preserved
mandel66_minor(arr) / array_from_mandel66(M)	6×6 Kelvin-Mandel ↔ 3×3×3×3, minor-symmetrizing, no major symmetry assumed
ti8_params_from_KM(M) / KM_from_ti8_params(p)	the eight Walpole coefficients of an axially-invariant tensor about e₃, read off exactly — the non-major-symmetric counterpart of TensND.ti_params_from_KM
ti_average_mandel66(M, n) / iso_average_mandel66(M)	the same averages, on a 6×6 block rather than on a tensor (the ageing-viscoelastic Volterra path)

See the TensND documentation for the full docstrings.