API — Interactions

MeanFieldHomogenization.Interactions Module

MeanFieldHomogenization.Interactions

Two-inclusion interaction tensors — the shared numerical ingredient of the N-body homogenization models of the package.

Where every one-site scheme (Dilute, Mori-Tanaka, Maxwell, PCW, self-consistent, differential) needs only the Hill tensor of a single inclusion in a reference medium, the equivalent inclusion method (EquivalentInclusion) and the cluster model (ClusterModel) need one object more: the tensor measuring (minus) the field induced in one inclusion by a uniform polarization of another, whose self term is   . Brisard et al. (2014) (2014), §3.1, observe that their influence pseudotensors of order zero coincide with the interaction tensors of Berveiller et al. (1987) and Molinari & El Mouden (1996) — so the two families share this module.

Contents

• api.jl : interaction_tensor, self_interaction_tensor and the _pair_kernel table (plus the sign convention, which is worth reading before transcribing formulas)

• pair_ball_iso.jl : closed forms for two balls (3D) or two disks (2D) in an isotropic reference — exact at any separation

• pair_multipole.jl : truncated multipole expansion for general ellipsoids

• pair_quadrature.jl : brute-force product quadrature, the validation oracle

• lattice_sums.jl : sums over the periodic images of a source inclusion

Two-inclusion interaction tensors

MeanFieldHomogenization.Interactions.interaction_tensor Function

interaction_tensor(incl_a, incl_b, r, P₀; method=:auto, kw...) -> AbstractTens

Interaction tensor between two non-overlapping inclusions embedded in an infinite reference medium P₀, with r the vector joining the center of incl_a to the center of incl_b.

It is the double average of the real-space Green operator over the two regions,

so that contracting it with a uniform polarization carried by incl_b returns minus the average field induced in incl_a — strain in elasticity, gradient of the temperature in conduction. Its self counterpart is self_interaction_tensor, equal to .

P₀ may be a 4th-order stiffness (elasticity) or a 2nd-order conductivity tensor; dispatch selects the corresponding formulation, in 2D or 3D. Two balls (3D) or two disks (2D) in an isotropic reference are evaluated by the closed form of Molinari and El Mouden (1996) and (Berveiller et al., 1987), which is exact at any separation; other geometries use the truncated multipole expansion of (Brisard et al., 2014), §4.2.

method selects the back-end explicitly: :analytical, :multipole, :quadrature, or :auto (default).

This is the shared numerical ingredient of both N-body models in the package, EquivalentInclusion and ClusterModel — Brisard et al. (2014), §3.1, note that their influence pseudotensors of order    coincide with the interaction tensors of Molinari & El Mouden.

Two sign conventions exist

The package follows Brisard et al. (2023), for which   . Molinari and El Mouden (1996) and Berveiller et al. (1987) use   , so a formula taken from them — their Appendix A table in particular — must be flipped before it is compared with anything here.

An anisotropic reference is supported in 3D elasticity and in conduction; only plane-strain elasticity with an anisotropic P₀ is missing, and it raises an ArgumentError naming the limitation rather than falling back to an isotropic kernel.

MeanFieldHomogenization.Interactions.self_interaction_tensor Function

self_interaction_tensor(incl, P₀; kw...) -> AbstractTens

Self term of the interaction family,   : minus the average field induced in an inclusion by its own uniform polarization, which is the Eshelby result    .

It is the Hill polarization tensor, so it inherits every back-end of hill_tensor — closed forms for isotropic and transversely isotropic references, the residue algorithm, the DECUHR and nested-QuadGK cubatures, in 2D and 3D, for elasticity and for conduction. Keyword arguments are forwarded to hill_tensor.

That the self term is plus the Hill tensor is the whole reason the package follows Brisard et al. (2023) rather than Molinari's opposite sign: it makes the N-body kernels and the one-site schemes share one object, and it is why the cluster model collapses onto Mori-Tanaka when the cluster is reduced to a single inclusion (Molinari and El Mouden, 1996; App. C).

Periodic images

MeanFieldHomogenization.Interactions.lattice_interaction_tensor Function

lattice_interaction_tensor(incl_a, incl_b, r, P₀, L, R_c; kw...) -> AbstractTens

Sum of interaction_tensor over every periodic image of the source inclusion lying within the cluster radius R_c of the receiver, for a cubic cell of side L:

When incl_a and incl_b are the same inclusion of the cell ( ), the null translation is skipped: the self term is not part of this sum, it is self_interaction_tensor.

Returns zero of the appropriate order when no image falls inside the cutoff, which is the physically correct answer (a cluster reduced to the receiver alone) and is what makes the cluster model degenerate exactly onto Mori-Tanaka.

MeanFieldHomogenization.Interactions.periodic_images Function

periodic_images(r, L, R_c; skip_self=false) -> Vector

Translations   ,  , whose norm is at most , for a cubic period . With skip_self = true the null translation is dropped, which is what a receiver needs when summing over its own family of images.

Enumerating a ball rather than a box is what makes the truncated lattice sum converge — see the file header and (Molinari and El Mouden, 1996), App. B.

Real-space Green operator

MeanFieldHomogenization.Core.green_operator Function

green_operator(P₀, x; kw...) -> AbstractTens

Real-space Green operator of an infinite medium of reference property P₀, evaluated at   — the regular kernel of the Lippmann-Schwinger equation, which the package writes in the sign convention of Brisard et al. (2023):

Dispatches on the symmetry class: an isotropic reference goes to the closed form green_operator_iso, anything else to green_operator_aniso; green_nodes sets the quadrature order of the latter (the name is distinct from the nodes of interaction_tensor's own quadrature back-end, so both can be given at once). Elasticity and conduction, 2D and 3D, are all covered; the one gap is plane-strain elasticity with an anisotropic reference, whose Green function needs the Stroh formalism rather than the Barnett line integral, and which raises an ArgumentError.

This is the entry point the two-inclusion interaction kernels call — through its SArray-returning twin _green_operator, which is the same dispatch without the Tens wrapper — so lifting the isotropic restriction here lifts it for interaction_tensor as well.

MeanFieldHomogenization.Core.green_operator_iso Function

green_operator_iso(C₀::TensISO{4,3}, x) -> Tens{4,3}

Real-space Green operator of an infinite isotropic elastic matrix C₀, evaluated at  :

where is the Kelvin Green function and the brackets denote symmetrization with respect to and . The leading minus is the convention of Brisard et al. (2023) — see the file header.

With  ,   and   , the second gradient of the Kelvin solution reads

A uniform polarization carried by a volume around the origin induces the strain    at — exact in the far field, and exact at any separation when the source region is a ball, because each component of is then averaged over a region where it is harmonic.

Returned as a Tens{4,3} with both minor symmetries and the major symmetry. Throws a DomainError at the origin. Type-generic (Float64, ForwardDiff.Dual, symbolic scalars).

See also green_gradient_iso, dipole_displacement_iso.

green_operator_iso(K₀::TensISO{2,3}, x) -> Tens{2,3}

Conduction counterpart: minus the Hessian of the scalar Green function of an infinite isotropic medium of conductivity ,

It is traceless away from the origin — the isotropic part of the interaction between two distinct inclusions vanishes identically in conduction, exactly as the deviatoric statement    does in elasticity.

The order-2 objects follow the    convention of the package (see the notation page): the field this kernel transports is the gradient and the source is the flux polarization , so every sign matches the elastic case symbol for symbol.

green_operator_iso(K₀::TensISO{2,2}, x) -> Tens{2,2}

Two-dimensional conduction counterpart, from   :

Multiplied by the area of a source disk this is literally the plane form given by Brisard et al. (2023),     , which is the sharpest available check that the package and the paper share one convention.

green_operator_iso(C₀::TensISO{4,2}, x) -> Tens{4,2}

Plane-strain elastic Green operator of an infinite isotropic matrix. The two-dimensional Kelvin solution is

whose second gradient gives

with   . is the plane-strain Poisson ratio of the reference medium. The result is symmetrized — and negated, per the convention of the file header — exactly as in 3D.

Anisotropic reference media

MeanFieldHomogenization.Core.green_function_aniso Function

green_function_aniso(C₀, x; nodes = 32) -> Tens{2,3}

Displacement Green function of an infinite medium of arbitrary anisotropic stiffness C₀, evaluated at   by the Barnett line integral

with   and  . The integrand is smooth and periodic, so a Gauss-Legendre rule converges fast: the Green function itself is at machine accuracy by 16 nodes.

Reduces to the Kelvin solution for an isotropic C₀, which is how the implementation is tested. Type-generic, including ForwardDiff.Dual — the Green operator is obtained by differentiating this very function, through the SArray-returning kernel _green_function_aniso rather than through this Tens wrapper.

MeanFieldHomogenization.Core.green_operator_aniso Function

green_operator_aniso(C₀, x; nodes = 32) -> Tens{4,3}

Real-space Green operator of an arbitrary anisotropic elastic medium,

with the leading minus of the (Brisard et al., 2023) convention, exactly as in the isotropic closed form.

The second gradient is taken with ForwardDiff through the line integral of green_function_aniso — exact, and reusing one verified expression rather than a second hand-derived one.

Two derivatives cost accuracy, so the operator needs a finer rule than the function — and how much finer depends on the anisotropy. For an isotropic stiffness the Christoffel inverse is a low-order trigonometric polynomial and the rule is exact by 16 nodes. For a genuinely anisotropic one (a cubic stiffness with its shear constant detuned by a third) the measured relative error against a 128-node reference is 1.4e-2 at 16 nodes, 8e-4 at 24, 1.7e-5 at 32 and 1.5e-12 at 64.

The default of 32 sits well below the truncation error of the multipole expansion that consumes it, at half the cost of 64 (the cost is linear in the node count: about 50 µs per interaction tensor at 32 nodes on a desktop machine). Raise it for a strongly anisotropic reference.

Considerably more expensive than the isotropic closed form (green_operator_iso) — tens of microseconds against a fraction of one — which is why the dispatcher green_operator prefers the latter whenever the reference is isotropic.

green_operator_aniso(K₀::AbstractTens{2,3}, x) -> Tens{2,3}

Conduction counterpart in three dimensions. Here the Green function is elementary,

so minus its Hessian is written in closed form: with    and   ,

It is K₀-traceless in the sense   , the anisotropic form of the vanishing isotropic part.

green_operator_aniso(K₀::AbstractTens{2,2}, x) -> Tens{2,2}

Two-dimensional conduction counterpart, from   :

MeanFieldHomogenization.Core.gauss_legendre_nodes Function

gauss_legendre_nodes(n, lo = -1, hi = 1) -> (x, w)

Gauss-Legendre nodes and weights on [lo, hi], from the Golub-Welsch eigenvalue problem of the Jacobi matrix. A few lines, no extra dependency, and plenty accurate at the orders used in the package.

The rule for a given n is computed once and memoized — see the note above the definition; this is a load-bearing optimization, not a micro-tuning.