MeanFieldHomogenization.jl

Appearance

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

The whole page as one script

Every block of this page, in order, is scripts/80_custom_inclusion_contract.jl. From a clone of the repository, julia --project=docs -e 'using Pkg; Pkg.instantiate()' once, then julia scripts/80_custom_inclusion_contract.jl, runs it end to end.

Before this page

The tutorial Hill polarization tensors in practice, whose Hill tensor is the first of the three entry gates, and the theory page Concentration and contribution tensors, which derives every tensor of the contract from that one.

MeanFieldHomogenization 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 leveled: you implement the lowest gate you can reach and the package derives everything above it.

GateYou supplyYou get
Athe Hill tensor the 8 concentration tensors, the 4 contribution tensors, every scheme
Bthe concentration tensor (both sides of it if the inclusion is heterogeneous)the derived concentrations, the contributions, every scheme
Cthe contribution tensors , 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.

julia
using MeanFieldHomogenization
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.

julia
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 = MeanFieldHomogenization.inclusion_basis(ell)
TensND.RotatedBasis{3, Float64}
→ basis: 3×3 Matrix{Float64}:
  0.730682  -0.644218  0.226026
  0.615445   0.764842  0.190379
 -0.29552    0.0       0.955336
→ dual basis: 3×3 Matrix{Float64}:
  0.730682  -0.644218  0.226026
  0.615445   0.764842  0.190379
 -0.29552    0.0       0.955336
→ covariant metric tensor: 3×3 TensND.Id2{3, Float64}:
 1.0   ⋅    ⋅ 
  ⋅   1.0   ⋅ 
  ⋅    ⋅   1.0
→ contravariant metric tensor: 3×3 TensND.Id2{3, Float64}:
 1.0   ⋅    ⋅ 
  ⋅   1.0   ⋅ 
  ⋅    ⋅   1.0

Gate A — the Hill tensor ​

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

julia
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 concentration tensor ​

For a morphology whose has no convenient closed form, or when an external solver returns directly.

The inclusion here is homogeneous, so the strain-side tensor alone is enough: the stress side follows from   . A heterogeneous inclusion has no such and must supply flux_gradient_concentration as well, as LayeredSphere does. Omitting it is silent: Dilute and MoriTanaka stay right, only the self-consistent schemes drift. As for gate A, each callback serves both physics.

julia
gate_B = CustomInclusion(
    semi; basis,
    gradient_gradient_concentration = (P₁, P₀; kw...) -> gradient_gradient_concentration(ell, P₁, P₀; kw...),
)
CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{gradient_gradient_concentration::Main.var"#7#8"}}((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, (gradient_gradient_concentration = Main.var"#7#8"(),))

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 , so Mori-Tanaka and the self-consistent schemes stay out of reach.

julia
gate_C = CustomInclusion(
    semi; basis,
    property_contribution = (P₁, P₀; kw...) -> property_contribution(ell, P₁, P₀; kw...),
    inverse_property_contribution = (P₁, P₀; kw...) -> inverse_property_contribution(ell, P₁, P₀; kw...),
)
CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{property_contribution::Main.var"#13#14", inverse_property_contribution::Main.var"#15#16"}}((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, (property_contribution = Main.var"#13#14"(), inverse_property_contribution = Main.var"#15#16"()))

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.

julia
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 concentration and contribution tensors are derived

✓ CustomInclusion{3, Float64, TensND.RotatedBasis{3, Float64}, @NamedTuple{property_contribution::Main.var"#13#14", inverse_property_contribution::Main.var"#15#16"}} satisfies the inclusion contract (physics = :elasticity, amount = :fraction)
  ✓ level 1, gate C: `property_contribution` / `inverse_property_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 ​

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

c1111(C) = canonical_components(C)[1, 1, 1, 1]
k11(K) = canonical_components(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  , and there is nothing to implement. For a flat object measured by a density the prefactor is geometric — for an elliptical crack of Budiansky density   — and that is exactly the seam density_factor fills.

julia
crack = EllipticCrack(1.0, 0.25)
flat = CustomInclusion(
    (1.0, 0.25, 0.0);
    basis = MeanFieldHomogenization.inclusion_basis(crack),
    density_factor = 4π / 3,
    inverse_property_contribution = (P₀; kw...) -> inverse_property_contribution(crack, P₀; kw...),
    property_contribution = (P₀; kw...) -> property_contribution(crack, P₀; kw...),
)

function cracked(geom, ε)
    r = RVE(; distribution_shape = Ellipsoid(1.0))
    add_phase!(r, :M, Ellipsoid(1.0), Dict(:C => C₀); fraction = :rest)
    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).

julia
function oriented(geom_at, nbins, f)
    r = RVE(; distribution_shape = Ellipsoid(1.0))
    add_phase!(r, :M, Ellipsoid(1.0), Dict(:C => C₀); fraction = :rest)
    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 = MeanFieldHomogenization.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(canonical_components(C_ref)[3, 3, 3, 3], digits = 10)
)
println(
    "  custom  : ", typeof(C_cus).name.name, "   C₃₃₃₃ = ",
    round(canonical_components(C_cus)[3, 3, 3, 3], digits = 10)
)
@printf "  ‖Δ‖ = %.3e\n" norm(C_cus - C_ref)

Orientation distribution (12 polar bins) + TISymmetrize
------------------------------------------------------------------------------
  native  : TensTI   C₃₃₃₃ = 49.5682276961
  custom  : TensTI   C₃₃₃₃ = 49.5682276961
  ‖Δ‖ = 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.

Where to go next ​

The callbacks above return closed-form tensors, but the contract only asks for tensors. An inclusion whose response is a neural network supplies them from a trained surrogate, through gates A and B, and examines what stays exact whatever the quality of the fit. Morphologies solved by finite elements enter through the same gates, as described in Finite-element inclusions.

This page was generated using Literate.jl.