API — Constitutive

MeanFieldHomogenization.Constitutive Module

MeanFieldHomogenization.Constitutive

MeanFieldHomogenization as a constitutive law at each Gauss point of a structural finite-element computation — the role an MFront behavior or an Abaqus UMAT plays.

Not to be confused with MeanFieldHomogenization.FiniteElements

FiniteElements points the other way: it uses finite elements inside the package, to solve the Eshelby problem of a single inclusion whose response no closed form covers. Here the finite-element code is the caller, and a whole microstructure plays the part of one material point.

The contract

A material bundles a microstructure, a scheme, and whatever internal state the model needs. An FE driver builds it once, allocates one state per quadrature point, and calls material_response in its element loop:

mat    = HomogenizedElastic(rve, MoriTanaka())
states = [[initial_state(mat) for _ in 1:nqp] for _ in 1:ncells]
cache  = MaterialCache()

r  = material_response(mat, ε, states[cell][qp], Δt; cache = cache)
σ  = to_tensors(stress(r))       # -> Tensors.SymmetricTensor{2,3}, global frame
ℂ  = to_tensors(tangent(r))      # -> Tensors.SymmetricTensor{4,3}, global frame
states[cell][qp] = state(r)

The contract is multi-gradient / multi-flux — the shape MGIS uses for MFront's generic behaviors — because the poroelastic case takes a strain and a pore pressure and returns a stress and a variation of fluid content. A purely mechanical law only ever sees the one-gradient specialization.

Contents

AbstractMFHMaterial, AbstractMaterialState	the contract
material_response, initial_state	the two required methods
MaterialResponse, stress, tangent, state	what a response carries
HomogenizedElastic	linear law from any scheme — the control case
MicrocrackedMaterial	crack families that open and close; exact piecewise-linear tangent
FracturedPoroelasticRock	the saturated fractured rock: two gradients, two fluxes, evolving permeability
MaterialCache, cache_stats	memoization on the microstructural configuration
to_tensors, from_tensors	the Tensors.jl bridge
plane_strain_response	driving a 3-D material from a 2-D element loop
check_material_interface	conformance checker

Two rules that are easy to break

• Frames. A Tensors.jl tensor has no basis; its components are global. A TensND tensor returns its components in its own basis, and a homogenized property whose RVE carries tilted inclusions comes back rotated. Always cross the boundary through to_tensors / from_tensors, never through get_array.

• State is immutable. material_response returns a new state instead of mutating its argument, so a rejected Newton iteration is recovered by keeping the old one.

The contract

MeanFieldHomogenization.Constitutive.AbstractMFHMaterial Type

AbstractMFHMaterial

Supertype of every Gauss-point constitutive law built on a homogenization model. A concrete material bundles the microstructure (an RVE or any AbstractHomogenizationCell), the scheme used to upscale it, and whatever is needed to evaluate it repeatedly and cheaply.

The contract

A material must implement:

• initial_state — the internal state at the start of the computation;

• material_response — the incremental response.

and may implement:

• gradient_names, flux_names, tangent_blocks — self-description, which check_material_interface uses and which a generic FE driver can introspect instead of hard-coding field names.

Materials with a single gradient and a single flux (every purely mechanical law) need only the three-argument form of material_response; the multi-field machinery then costs them nothing.

See also HomogenizedElastic.

MeanFieldHomogenization.Constitutive.AbstractMaterialState Type

AbstractMaterialState

Supertype of the internal state carried by a Gauss-point material between steps: crack aperture ratios, open/closed flags, fracture conductivities, and so on.

A state is immutable by convention. material_response returns a new one rather than mutating its argument, which is what makes a rejected Newton iteration trivially recoverable — the caller simply keeps the old state — and what keeps the law usable from a multithreaded element loop. The FE driver holds two arrays, states and states_old, and swaps them once a step has converged.

NoState is the state of a law that has none.

MeanFieldHomogenization.Constitutive.NoState Type

NoState() <: AbstractMaterialState

The internal state of a law that carries none — a linear elastic material, for instance. Costs nothing to store per quadrature point.

MeanFieldHomogenization.Constitutive.MaterialResponse Type

MaterialResponse(fluxes, tangents, state)

What material_response returns.

• fluxes — a NamedTuple of the thermodynamic forces, e.g. (; σ) for a mechanical law, (; σ, φ) for a poroelastic one (φ being the variation of fluid content).

• tangents — a NamedTuple of the tangent blocks, keyed by flux then gradient: (; σε = ℂ, σp = -𝐁, φε = 𝐁, φp = 1/M). Keys concatenate the flux and gradient names so the block structure stays readable at the call site.

• state — the updated AbstractMaterialState.

Access the common mechanical case with stress and tangent rather than by digging into the fields.

MeanFieldHomogenization.Constitutive.initial_state Function

initial_state(m::AbstractMFHMaterial) -> AbstractMaterialState

The internal state a fresh quadrature point starts from. Called once per quadrature point when the FE driver allocates its state arrays.

MeanFieldHomogenization.Constitutive.material_response Function

material_response(m, gradients, state_old, Δt; cache = nothing) -> MaterialResponse
material_response(m, ε, state_old, Δt; cache = nothing) -> MaterialResponse

Integrate the constitutive law over one step and return the fluxes, the tangent blocks and the new state.

gradients is a NamedTuple — (; ε) for a mechanical law, (; ε, p) for a poroelastic one. The second form accepts a bare strain tensor and wraps it, so a mechanical law reads

r = material_response(mat, ε, st, Δt)
σ, ℂ, st_new = stress(r), tangent(r), state(r)

state_old is not mutated: the response carries a new state. That is what makes an FE Newton iteration that has to be retried safe, and what allows a threaded element loop.

Δt is the time increment; rate-independent laws ignore it. cache is an optional MaterialCache shared across quadrature points — for a homogenization-backed law it is the difference between a scheme solve per Gauss point and one per distinct microstructural state.

Tensors are TensND objects in the global frame

ε is expected in the canonical (global) frame, which is where an FE code computes it. Use from_tensors on the way in and to_tensors on the way out; do not hand raw component arrays across the boundary.

MeanFieldHomogenization.Constitutive.stress Function

stress(r::MaterialResponse)

The σ flux of a response — the stress of a mechanical or poroelastic law.

MeanFieldHomogenization.Constitutive.tangent Function

tangent(r::MaterialResponse)

The σε tangent block — the consistent tangent stiffness an FE code assembles into its Jacobian.

MeanFieldHomogenization.Constitutive.state Function

state(r::MaterialResponse)

The updated internal state carried by a response.

Self-description

MeanFieldHomogenization.Constitutive.gradient_names Function

gradient_names(m::AbstractMFHMaterial) -> Tuple{Vararg{Symbol}}

The gradients the law consumes, in order. Defaults to (:ε,).

MeanFieldHomogenization.Constitutive.flux_names Function

flux_names(m::AbstractMFHMaterial) -> Tuple{Vararg{Symbol}}

The thermodynamic forces the law produces, in order. Defaults to (:σ,).

MeanFieldHomogenization.Constitutive.tangent_blocks Function

tangent_blocks(m::AbstractMFHMaterial) -> Tuple{Vararg{Symbol}}

The tangent blocks the law supplies, keyed flux then gradient. Defaults to (:σε,).

A generic FE driver can read this to know which blocks to assemble instead of hard-coding them, which is what makes the same driver serve a mechanical and a poroelastic material.

MeanFieldHomogenization.Constitutive.transport_property Function

transport_property(m::AbstractMFHMaterial, state) -> Union{Nothing,AbstractTens{2,3}}

The 2nd-order transport property (permeability, conductivity, diffusivity) implied by the current internal state, or nothing for a material that carries none.

This is not a flux: it is a coefficient the surrounding balance equation needs — Darcy's law, in the fractured-reservoir case, where the permeability follows the fracture apertures and so changes with the state. Keeping it off the flux/tangent structure reflects that it belongs to a different balance equation than the one the material closes.

transport_property(m::FracturedPoroelasticRock, st) -> Union{Nothing,Tens{2,3}}

Effective permeability implied by the current apertures, through fracture_permeability.

nothing when no family is a ConductiveCrack: the material is then purely poroelastic and the flow problem is not its business.

MeanFieldHomogenization.Constitutive.check_material_interface Function

check_material_interface(m::AbstractMFHMaterial; verbose = true, kw...) -> Bool

Exercise a Gauss-point material against the contract and report what it does and does not honor. Returns true when every check passes.

What is verified:

1. initial_state returns an AbstractMaterialState;

2. material_response accepts both the bare-strain and the NamedTuple forms and returns a MaterialResponse;

3. every name declared by flux_names and tangent_blocks is actually present in the response;

4. the response is frame-honest: the fluxes and tangents come back in the canonical basis, which is what an FE code will read them as;

5. state_old is not mutated — the contract requires a new state, so that a rejected Newton iteration can simply be retried;

6. the declared σε block agrees with a central finite difference of the stress. This is the check that matters: a plausible-but-wrong tangent costs quadratic convergence and nothing else, so it is easy to ship.

The tangent check is skipped where it is meaningless — at a state on a non-smooth branch of the law. Pass probe to move the probe strain somewhere smooth.

Keywords

• probe — the strain at which to exercise the law (default: a small non-symmetric-looking triaxial state, so a law that only works in uniaxial tension is caught).

• Δ — finite-difference step (default 1e-7).

• rtol — tolerance on the tangent comparison (default 1e-5).

• verbose — print the report (default true).

julia> check_material_interface(HomogenizedElastic(rve, MoriTanaka()));

Materials

MeanFieldHomogenization.Constitutive.HomogenizedElastic Type

HomogenizedElastic(cell, scheme; property = :C, kw...)
HomogenizedElastic(C_hom)

Linear elastic Gauss-point material whose stiffness comes from a homogenization scheme.

The scheme is run once, at construction, and the resulting C_hom is stored in the canonical frame — an RVE with oriented inclusions returns its estimate in a rotated basis, and a material handed to an FE code must speak the global frame (see to_tensors).

rve = RVE()
add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => TensISO{3}(3 * 30.0, 2 * 18.0)); fraction = :rest)
add_phase!(rve, :I, Ellipsoid(1.0), Dict(:C => TensISO{3}(3 * 90.0, 2 * 60.0)); fraction = 0.2)

mat = HomogenizedElastic(rve, MoriTanaka())
r   = material_response(mat, ε, initial_state(mat), 0.0)
σ, ℂ = stress(r), tangent(r)

The second form takes a stiffness directly, which is useful for a control run or when the upscaling was done elsewhere.

Extra keyword arguments are forwarded to homogenize, so solver tolerances travel with the material: HomogenizedElastic(rve, SelfConsistent(); abstol = 1e-12).

MeanFieldHomogenization.Constitutive.stiffness Function

stiffness(m::HomogenizedElastic) -> Tens{4,3}

The homogenized stiffness the material was built with, in the canonical frame.

MeanFieldHomogenization.Constitutive.MicrocrackedMaterial Type

MicrocrackedMaterial(rve, scheme; families, ω₀, kw...)

Gauss-point material for a solid holding crack families that open and close with the loading.

• rve — the microstructure, whose crack phases are the families. Their densities, normals and radii are fixed data; only the apertures evolve.

• scheme — any scheme supporting the per-family decomposition (SelfConsistent, MoriTanaka).

• families — the crack phase names, in order. Defaults to every CrackDensity phase of rve.

• ω₀ — initial aspect ratio of each family. A family closes when its aspect ratio reaches zero, so ω₀ sets how much compression it tolerates.

The law is piecewise linear: within a branch the tangent is exactly of the current open set, and a step crossing a closure is split at the crossing.

rve = RVE()
add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => C₀); fraction = :rest)
add_phase!(rve, :F1, PennyCrack(1.0), Dict(:C => C₀); density = 0.1)

mat = MicrocrackedMaterial(rve, MoriTanaka(); ω₀ = (1.0e-3,))
r   = material_response(mat, ε, initial_state(mat), 0.0; cache = MaterialCache())

Closed means invisible, not welded

A closed family is dropped from the homogenization, i.e. it transmits both normal and tangential tractions. Frictional sliding on closed faces is a different model and is not what this material does.

See also HomogenizedElastic for the linear control case, and crack_family_compliances for the decomposition the aperture update rests on.

MeanFieldHomogenization.Constitutive.CrackedState Type

CrackedState(ω, open, ε, σ)

Internal state of a MicrocrackedMaterial: the current aspect ratio ω[i] of each crack family, whether it is open[i], and the strain ε and stress σ the state was reached at.

The law is incremental, so both ε and σ are carried: ε to form the strain increment that drives the apertures, and σ because the stress must be integrated along the path rather than recomputed as  . Those two differ as soon as a family closes — the medium stiffens from that point on, and only from that point on, so a total-strain formula would make the stress jump discontinuously at closure.

MeanFieldHomogenization.Constitutive.open_set Function

open_set(st::CrackedState) -> NTuple{N,Bool}

Which crack families are currently open. This tuple — and nothing else about the state — is what the homogenized stiffness depends on, which is what makes MaterialCache effective: at most 2^N distinct stiffnesses exist, however many quadrature points there are.

MeanFieldHomogenization.Constitutive.apertures Function

apertures(st::CrackedState) -> NTuple{N,T}

Current aspect ratios   of the crack families.

MeanFieldHomogenization.Constitutive.FracturedPoroelasticRock Type

FracturedPoroelasticRock(rve, scheme; ω₀, C₀, k_matrix, porosity_ref, kw...)

The saturated fractured rock of [71]: a Gauss-point material with two gradients and two fluxes , whose fractures open and close and whose permeability follows their apertures.

• rve — matrix plus crack families. Families given as ConductiveCrack also carry the hydraulic side; ordinary cracks make the material purely poroelastic and transport_property returns nothing.

• ω₀ — initial aspect ratios, as for MicrocrackedMaterial.

• k_matrix — matrix conductivity, small but non-zero (see fracture_permeability).

The tangent blocks are :σε , :σp  , :φε and :φp , all recomputed whenever the open/closed set changes and cached on it.

mat = FracturedPoroelasticRock(rve, SelfConsistent(); ω₀ = (1.0e-4,), k_matrix = 1.0e-18)
r = material_response(mat, (; ε = ε, p = p), st, Δt; cache = cache)
r.fluxes.σ, r.fluxes.φ, r.tangents.σp, transport_property(mat, r.state)

Incompressible fluid

comes from inverse_biot_modulus, which assumes   — the setting of the paper. See that docstring for the compressible generalization.

MeanFieldHomogenization.Constitutive.PoroFracturedState Type

PoroFracturedState(mech, C)

Internal state of a FracturedPoroelasticRock: the mechanical state of the fracture families (CrackedState) and their current conductivities C[i].

p is the pore pressure the state was reached at, kept for the same reason the mechanical state keeps its strain: the driver is an increment.

The conductivities are carried separately because they follow the apertures by the cubic law   and feed a different balance equation — Darcy's, through transport_property — rather than the mechanical one.

MeanFieldHomogenization.Constitutive.conductivities Function

Current fracture conductivities of a PoroFracturedState.

MeanFieldHomogenization.Constitutive.fluid_content Function

fluid_content(m::FracturedPoroelasticRock, st; cache = nothing) -> Real

Fluid-content variation     at the state st, from the natural initial state.

A transient flow problem discretizes , so its residual needs at both ends of the step. material_response returns the new one; this returns the old one, recomputed from the stored strain and pressure with the poroelastic parameters of that state's own open set — which is why a driver never has to carry alongside the state.

Δφ = material_response(m, (; ε, p), st_old, Δt).fluxes.φ - fluid_content(m, st_old)

Memoization

MeanFieldHomogenization.Constitutive.MaterialCache Type

MaterialCache()

Memoization store shared by the quadrature points of a AbstractMFHMaterial.

Pass one to material_response via its cache keyword and reuse it across the whole element loop — the entries are keyed on the microstructural configuration, not on the point, so every quadrature point sharing a configuration pays for it once.

cache = MaterialCache()
for cell in CellIterator(dh), qp in 1:getnquadpoints(cv)
    r = material_response(mat, ε[qp], states[qp], Δt; cache = cache)
end
@show cache_stats(cache)     # (; hits, misses, entries)

Not thread-safe

A MaterialCache is a plain Dict with no lock. Use one cache per thread in a threaded element loop, or none at all. Sharing one across threads is a data race, and the failure mode — a torn read of a partially built entry — would surface as a wrong stiffness rather than as an error.

See also cached!, cache_stats.

MeanFieldHomogenization.Constitutive.cached! Function

cached!(f, cache, key)
cached!(f, ::Nothing, key)

Return cache.entries[key], computing it with f() on a miss. With nothing in place of a cache, f() is always called — so a material can be written once and used with or without memoization.

MeanFieldHomogenization.Constitutive.cache_stats Function

cache_stats(cache) -> NamedTuple

(; hits, misses, entries) for a MaterialCache, or a zeroed tuple for nothing.

Worth looking at once on a real mesh: for a model whose fracture network is uniform over a region, entries should settle at the number of distinct open/closed configurations actually visited (a handful), not grow with the number of quadrature points. If it grows with the mesh, the cache key is carrying continuous data and the memoization is not doing its job.

MeanFieldHomogenization.Constitutive.reset_cache! Function

reset_cache!(cache) -> MaterialCache

Empty a MaterialCache and zero its counters. Needed whenever the underlying microstructure changes in a way the key does not capture — a different fracture network per material region, for instance, if the network id is not part of the key.

Tensors.jl bridge

MeanFieldHomogenization.Constitutive.to_tensors Function

to_tensors(t::AbstractTens{2,3}) -> Tensors.SymmetricTensor{2,3}
to_tensors(t::AbstractTens{4,3}) -> Tensors.SymmetricTensor{4,3}

Convert a TensND tensor to its Tensors.jl counterpart, expressed in the global (canonical) frame.

This is the direction an FE code consumes: a stress, or a tangent stiffness, ready to be contracted with shape-function gradients.

Never bypass this with get_array

get_array(t) and t[i,j,k,l] return the components of t in its own basis. Any homogenized property whose RVE carries oriented inclusions — tilted crack families, in particular — is returned in a rotated basis, so those components are not the global ones. to_tensors performs the change of basis; reading the array directly silently rotates the material.

The 4th-order conversion assumes minor symmetry, which every stiffness, compliance and localization tensor in the package satisfies. Major symmetry is not assumed: Tensors.SymmetricTensor{4,3} carries 36 independent components and can hold a concentration tensor unchanged.

MeanFieldHomogenization.Constitutive.from_tensors Function

from_tensors(t::Tensors.SymmetricTensor{2,3}) -> Tens{2,3}
from_tensors(t::Tensors.SymmetricTensor{4,3}) -> Tens{4,3}

Convert a Tensors.jl tensor coming from an FE code into a TensND tensor in the canonical basis — which is what its components mean.

This is the direction an FE code produces: the strain at a quadrature point, assembled from shape-function gradients in the global frame.

MeanFieldHomogenization.Constitutive.plane_strain_response Function

plane_strain_response(m, ε₂, state_old, Δt; cache = nothing)

Drive a three-dimensional AbstractMFHMaterial from a plane-strain two-dimensional element loop.

ε₂ is the in-plane strain as a Tensors.SymmetricTensor{2,2} — what a 2-D Ferrite/Gridap element produces. Returns a NamedTuple:

field
σ	in-plane stress, SymmetricTensor{2,2}
C	in-plane tangent, SymmetricTensor{4,2}
σ₃₃	out-of-plane stress — not zero in plane strain
state	the updated internal state

everything in the global frame, ready to assemble.

The reduction is exact: plane strain means    , so the 3-D strain is the in-plane one padded with zeros, and the 2-D tangent is the in-plane block   of the 3-D one.

r = plane_strain_response(mat, ε₂, states[c][q], Δt; cache = cache)
Ke .+= ... r.C ...
states[c][q] = r.state

Plane strain only — not plane stress

Plane stress ( ) requires condensing the out-of-plane strain out of the law, which for a general anisotropic couples all six components and, for a material with internal state, has to be solved at every quadrature point. Reusing this function for plane stress would silently impose   instead of  , which is a different problem — note the σ₃₃ this returns is generally non-zero.

Anisotropy is not checked

A microstructure whose axes are not aligned with the plane produces a coupling in-plane and out-of-plane components ( and friends). Plane strain remains exact — those couplings only feed σ₃₃ and the out-of-plane shears, which the 2-D momentum balance does not see — but the resulting plane problem is then not the one a 2-D intuition expects.

MeanFieldHomogenization.Constitutive.voigt_stress Function

voigt_stress(σ) -> NTuple{6}
voigt_strain(ε) -> NTuple{6}

Components of a symmetric 2nd-order tensor in the Voigt convention used by Abaqus-style interfaces, ordered (11, 22, 33, 12, 13, 23), in the global frame.

The two differ by the engineering-shear factor: voigt_strain doubles the off-diagonal terms (γ₁₂ = 2ε₁₂), voigt_stress does not. Keeping them as separate functions is deliberate — a single voigt helper is the classic way to lose a factor of two between the strain that goes in and the stress that comes out.

Note this is not the Kelvin-Mandel convention used internally by TensND and Tensors.jl (which carries √2 instead, and is an isometry); Voigt appears here only at the boundary with codes that demand it, Abaqus-style UMATs above all.

MeanFieldHomogenization.Constitutive.voigt_strain Function

voigt_stress(σ) -> NTuple{6}
voigt_strain(ε) -> NTuple{6}

Components of a symmetric 2nd-order tensor in the Voigt convention used by Abaqus-style interfaces, ordered (11, 22, 33, 12, 13, 23), in the global frame.

The two differ by the engineering-shear factor: voigt_strain doubles the off-diagonal terms (γ₁₂ = 2ε₁₂), voigt_stress does not. Keeping them as separate functions is deliberate — a single voigt helper is the classic way to lose a factor of two between the strain that goes in and the stress that comes out.

Note this is not the Kelvin-Mandel convention used internally by TensND and Tensors.jl (which carries √2 instead, and is an isometry); Voigt appears here only at the boundary with codes that demand it, Abaqus-style UMATs above all.

MeanFieldHomogenization.Constitutive.stress_from_voigt Function

stress_from_voigt(v) -> Tens{2,3}

Inverse of voigt_stress: rebuild a stress tensor from (σ₁₁, σ₂₂, σ₃₃, σ₁₂, σ₁₃, σ₂₃), with no shear factor.

MeanFieldHomogenization.Constitutive.strain_from_voigt Function

strain_from_voigt(v) -> Tens{2,3}

Inverse of voigt_strain: rebuild a strain tensor from the six Voigt components (ε₁₁, ε₂₂, ε₃₃, γ₁₂, γ₁₃, γ₂₃), halving the engineering shears.