GitHub

The custom-inclusion contract: three entry gates, one answer

MeanFieldHom lets an arbitrary morphology take part in every homogenization scheme, in elasticity and in transport, without touching the package. This is the counterpart of the user_inclusion mechanism of the C++/Python echoes codebase, but levelled: you implement the lowest gate you can reach and the package derives everything above it.

GateYou supplyYou get
Athe Hill tensor $\mathbb P$the 8 localization tensors, the 4 contribution tensors, every scheme
Bthe localization tensor $\mathbb A$ (both sides of it if the inclusion is heterogeneous)the derived localizations, the contributions, every scheme
Cthe contribution tensors $\mathbb N$, $\mathbb H$the contribution-based schemes

This script drives the same physical inclusion through all three gates and checks that the schemes cannot tell them apart. The point is not the numbers — they are identical by construction — but that the plumbing is identical.

See the manual page Custom inclusions and the developer page Adding a new inclusion for the full contract.

using MeanFieldHom
using TensND
using Printf
using LinearAlgebra

The reference: a tilted prolate spheroid

Tilted, so that the effective tensor is genuinely anisotropic and a dropped rotation cannot hide behind an isotropic answer.

C₀ = TensISO{3}(3 * 30.0, 2 * 10.0)      # matrix
C₁ = TensISO{3}(3 * 60.0, 2 * 20.0)      # inclusion
K₀ = TensISO{3}(2.0)                     # matrix conductivity
K₁ = TensISO{3}(7.0)

ell = Ellipsoid(3.0, 1.0, 1.0; euler_angles = (0.3, 0.7, 0.0))
semi = (3.0, 1.0, 1.0)
basis = MeanFieldHom.inclusion_basis(ell)
3×3 TensND.RotatedBasis{3, Float64}:
  0.730682  -0.644218  0.226026
  0.615445   0.764842  0.190379
 -0.29552    0.0       0.955336

Gate A — the Hill tensor

One callback. Everything else is derived: the eight localization tensors, the four contribution tensors, and every scheme.

gate_A = CustomInclusion(
    semi; basis,
    hill_tensor = (P₀; kw...) -> hill_tensor(ell, P₀; kw...),
)
CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{hill_tensor::Main.var"#3#4"}}((3.0, 1.0, 1.0), [0.7306816499355124 -0.644217687237691 0.22602632124962302; 0.6154446635582734 0.7648421872844885 0.19037934406737264; -0.29552020666133955 0.0 0.955336489125606], true, nothing, (hill_tensor = Main.var"#3#4"(),))

P₀ is a 4th-order stiffness in elasticity and a 2nd-order conductivity in transport, so a single callback serves both physics.

Gate B — the localization tensor

For a morphology whose $\mathbb P$ has no convenient closed form, or when an external solver returns $\mathbb A$ directly.

The inclusion here is homogeneous, so the strain-side tensor alone is enough: the stress side follows from $\mathbb A_{\sigma\varepsilon} = \mathbb C_1:\mathbb A_{\varepsilon\varepsilon}$. A heterogeneous inclusion has no such $\mathbb C_1$ and must supply stress_strain_loc (and flux_gradient_loc) as well — LayeredSphere implements all four. Omitting them is silent: Dilute and MoriTanaka stay right, only the self-consistent schemes drift.

gate_B = CustomInclusion(
    semi; basis,
    strain_strain_loc = (P₁, P₀; kw...) -> strain_strain_loc(ell, P₁, P₀; kw...),
    gradient_gradient_loc = (P₁, P₀; kw...) -> gradient_gradient_loc(ell, P₁, P₀; kw...),
)
CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{strain_strain_loc::Main.var"#9#10", gradient_gradient_loc::Main.var"#11#12"}}((3.0, 1.0, 1.0), [0.7306816499355124 -0.644217687237691 0.22602632124962302; 0.6154446635582734 0.7648421872844885 0.19037934406737264; -0.29552020666133955 0.0 0.955336489125606], true, nothing, (strain_strain_loc = Main.var"#9#10"(), gradient_gradient_loc = Main.var"#11#12"()))

Gate C — the contribution tensors

The last resort, and the one that buys the least: contribution tensors alone do not determine the dilute concentration tensor $\mathbb A$, so Mori-Tanaka and the self-consistent schemes stay out of reach.

gate_C = CustomInclusion(
    semi; basis,
    stiffness_contribution = (P₁, P₀; kw...) -> stiffness_contribution(ell, P₁, P₀; kw...),
    compliance_contribution = (P₁, P₀; kw...) -> compliance_contribution(ell, P₁, P₀; kw...),
    conductivity_contribution = (P₁, P₀; kw...) -> conductivity_contribution(ell, P₁, P₀; kw...),
    resistivity_contribution = (P₁, P₀; kw...) -> resistivity_contribution(ell, P₁, P₀; kw...),
)
CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{stiffness_contribution::Main.var"#21#22", compliance_contribution::Main.var"#23#24", conductivity_contribution::Main.var"#25#26", resistivity_contribution::Main.var"#27#28"}}((3.0, 1.0, 1.0), [0.7306816499355124 -0.644217687237691 0.22602632124962302; 0.6154446635582734 0.7648421872844885 0.19037934406737264; -0.29552020666133955 0.0 0.955336489125606], true, nothing, (stiffness_contribution = Main.var"#21#22"(), compliance_contribution = Main.var"#23#24"(), conductivity_contribution = Main.var"#25#26"(), resistivity_contribution = Main.var"#27#28"()))

Checking conformance before wiring anything

check_inclusion_interface reports the level-0 methods, the gate it found, and — this is the useful part — what that gate leaves unavailable.

println("="^78)
check_inclusion_interface(gate_A)
println()
check_inclusion_interface(gate_C)
println("="^78)
==============================================================================
✓ CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{hill_tensor::Main.var"#3#4"}} satisfies the inclusion contract (physics = :elasticity, amount = :fraction)
  ✓ level 1, gate A: `hill_tensor` — all localization and contribution tensors are derived

✓ CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{stiffness_contribution::Main.var"#21#22", compliance_contribution::Main.var"#23#24", conductivity_contribution::Main.var"#25#26", resistivity_contribution::Main.var"#27#28"}} satisfies the inclusion contract (physics = :elasticity, amount = :fraction)
  ✓ level 1, gate C: `stiffness_contribution` / `compliance_contribution` supplied directly
    ! gate C alone does not provide the dilute concentration tensor, so `MoriTanaka`, `SelfConsistent` and `AsymmetricSelfConsistent` are unavailable for this volume-fraction phase (Voigt, Reuss, Dilute, DiluteDual, Maxwell, PonteCastanedaWillis and DifferentialScheme work).
==============================================================================

The schemes cannot tell them apart

function rve_with(geom, prop_m, prop_i, key)
    r = RVE(:M)
    add_matrix!(r, Ellipsoid(1.0), Dict(key => prop_m))
    add_phase!(r, :I, geom, Dict(key => prop_i); fraction = 0.25)
    return r
end

c1111(C) = get_array(C)[1, 1, 1, 1]
k11(K) = get_array(K)[1, 1]

ALL = (
    "Dilute" => Dilute(), "DiluteDual" => DiluteDual(),
    "MoriTanaka" => MoriTanaka(), "Maxwell" => Maxwell(),
    "PCW" => PonteCastanedaWillis(), "SelfConsistent" => SelfConsistent(),
    "ASC" => AsymmetricSelfConsistent(), "Differential" => DifferentialScheme(),
)
GATE_C_OK = ("Dilute", "DiluteDual", "Maxwell", "PCW", "Differential")

println("\nElasticity — C_eff[1,1,1,1]\n", "-"^78)
@printf "%-16s %14s %14s %14s %14s\n" "scheme" "Ellipsoid" "gate A" "gate B" "gate C"
for (name, sch) in ALL
    ref = c1111(homogenize(rve_with(ell, C₀, C₁, :C), sch, :C))
    vA = c1111(homogenize(rve_with(gate_A, C₀, C₁, :C), sch, :C))
    vB = c1111(homogenize(rve_with(gate_B, C₀, C₁, :C), sch, :C))
    if name in GATE_C_OK
        vC = c1111(homogenize(rve_with(gate_C, C₀, C₁, :C), sch, :C))
        @printf "%-16s %14.8f %14.8f %14.8f %14.8f\n" name ref vA vB vC
    else
        @printf "%-16s %14.8f %14.8f %14.8f %14s\n" name ref vA vB "n/a"
    end
end

println("\nConduction — K_eff[1,1]\n", "-"^78)
@printf "%-16s %14s %14s %14s %14s\n" "scheme" "Ellipsoid" "gate A" "gate B" "gate C"
for (name, sch) in (
        "Dilute" => Dilute(), "DiluteDual" => DiluteDual(),
        "MoriTanaka" => MoriTanaka(), "Maxwell" => Maxwell(),
    )
    ref = k11(homogenize(rve_with(ell, K₀, K₁, :K), sch, :K))
    vA = k11(homogenize(rve_with(gate_A, K₀, K₁, :K), sch, :K))
    vB = k11(homogenize(rve_with(gate_B, K₀, K₁, :K), sch, :K))
    if name in GATE_C_OK
        vC = k11(homogenize(rve_with(gate_C, K₀, K₁, :K), sch, :K))
        @printf "%-16s %14.8f %14.8f %14.8f %14.8f\n" name ref vA vB vC
    else
        @printf "%-16s %14.8f %14.8f %14.8f %14s\n" name ref vA vB "n/a"
    end
end

Elasticity — C_eff[1,1,1,1]
------------------------------------------------------------------------------
scheme                Ellipsoid         gate A         gate B         gate C
Dilute              50.34089121    50.34089121    50.34089121    50.34089121
DiluteDual          51.73523441    51.73523441    51.73523441    51.73523441
MoriTanaka          50.99598143    50.99598143    50.99598143            n/a
Maxwell             51.11639531    51.11639531    51.11639531    51.11639531
PCW                 51.11639531    51.11639531    51.11639531    51.11639531
SelfConsistent      51.07400466    51.07400466    51.07400466            n/a
ASC                 51.14371392    51.14371392    51.14371392            n/a
Differential        51.06857736    51.06857736    51.06857736    51.06857736

Conduction — K_eff[1,1]
------------------------------------------------------------------------------
scheme                Ellipsoid         gate A         gate B         gate C
Dilute               2.80034593     2.80034593     2.80034593     2.80034593
DiluteDual           3.42310410     3.42310410     3.42310410     3.42310410
MoriTanaka           2.87178043     2.87178043     2.87178043            n/a
Maxwell              2.93327474     2.93327474     2.93327474     2.93327474

Every column agrees to the last digit: the gates are algebraically equivalent reformulations of the same information, and the schemes consume whichever is available.

Level 2 — the amount × contribution seam

A scheme adds up an amount times a size-independent contribution. For a volume fraction that is simply $f\,\mathbb N$, and there is nothing to implement. For a flat object measured by a density the prefactor is geometric — $4\pi/3$ for an elliptical crack of Budiansky density $\varepsilon^{3\mathrm d} = N a b^2$ — and that is exactly the seam density_factor fills.

crack = EllipticCrack(1.0, 0.25)
flat = CustomInclusion(
    (1.0, 0.25, 0.0);
    basis = MeanFieldHom.inclusion_basis(crack),
    density_factor = 4π / 3,
    compliance_contribution = (P₀; kw...) -> compliance_contribution(crack, P₀; kw...),
    stiffness_contribution = (P₀; kw...) -> stiffness_contribution(crack, P₀; kw...),
)

function cracked(geom, ε)
    r = RVE(:M)
    add_matrix!(r, Ellipsoid(1.0), Dict(:C => C₀))
    add_phase!(r, :cr, geom, Dict(:C => C₀); density = ε)
    return r
end

println("\nFlat object with a density amount — C_eff[1,1,1,1]\n", "-"^78)
@printf "%-8s %18s %18s\n" "ε" "EllipticCrack" "CustomInclusion"
for ε in (0.0, 0.02, 0.05, 0.1)
    @printf "%-8.2f %18.10f %18.10f\n" ε c1111(homogenize(cracked(crack, ε), MoriTanaka(), :C)) c1111(homogenize(cracked(flat, ε), MoriTanaka(), :C))
end

Flat object with a density amount — C_eff[1,1,1,1]
------------------------------------------------------------------------------
ε             EllipticCrack    CustomInclusion
0.00          43.3333333333      43.3333333333
0.02          41.0671881033      41.0671881033
0.05          38.8743485537      38.8743485537
0.10          36.7513143409      36.7513143409

What is free: orientation averaging

IsoSymmetrize / TISymmetrize are applied by the scheme, after the callback returns — so a custom inclusion inherits them at no cost. The only obligations are to return tensors in the global frame, and to know that under symmetrization the scheme hands the kernel a pre-projected reference medium (isotropic by default).

function oriented(geom_at, nbins, f)
    r = RVE(:M)
    add_matrix!(r, Ellipsoid(1.0), Dict(:C => C₀))
    for (i, bin) in enumerate(polar_orientation_bins(nbins))
        add_phase!(
            r, Symbol(:I, i), geom_at(bin.θ), Dict(:C => C₁);
            fraction = f * bin.weight, symmetrize = TISymmetrize((0.0, 0.0, 1.0))
        )
    end
    return r
end

spheroid_at(θ) = Ellipsoid(3.0, 1.0, 1.0; euler_angles = (θ, 0.0, 0.0))
custom_at(θ) = let e = spheroid_at(θ)
    CustomInclusion(
        semi; basis = MeanFieldHom.inclusion_basis(e),
        hill_tensor = (P₀; kw...) -> hill_tensor(e, P₀; kw...)
    )
end

C_ref = homogenize(oriented(spheroid_at, 12, 0.2), MoriTanaka(), :C)
C_cus = homogenize(oriented(custom_at, 12, 0.2), MoriTanaka(), :C)

println("\nOrientation distribution (12 polar bins) + TISymmetrize\n", "-"^78)
println(
    "  native  : ", typeof(C_ref).name.name, "   C₃₃₃₃ = ",
    round(get_array(C_ref)[3, 3, 3, 3], digits = 10)
)
println(
    "  custom  : ", typeof(C_cus).name.name, "   C₃₃₃₃ = ",
    round(get_array(C_cus)[3, 3, 3, 3], digits = 10)
)
@printf "  max |Δ| = %.3e\n" maximum(abs, get_array(C_cus) .- get_array(C_ref))

Orientation distribution (12 polar bins) + TISymmetrize
------------------------------------------------------------------------------
  native  : TensTI   C₃₃₃₃ = 49.5682276961
  custom  : TensTI   C₃₃₃₃ = 49.5682276961
  max |Δ| = 0.000e+00

The exact transverse-isotropic average is applied identically in both cases: the custom inclusion never had to know that orientation averaging existed.

This page was generated using Literate.jl.