Skip to content

Ferrite backend

The material contract knows nothing about any finite-element library, so coupling to one needs almost no adapter. What Ferrite does not provide is per-quadrature-point material state, so that bookkeeping — and an element routine — is all this extension adds.

Activated by import Ferrite alone:

julia
import Ferrite
const EXT = Base.get_extension(
    MeanFieldHomogenization, :MeanFieldHomogenizationFerriteMaterialExt
)

Why a second extension

MeanFieldHomogenizationFerriteExt serves the opposite coupling and needs FerriteGmsh and Gmsh too. A structural computation that only wants a homogenized material law has no reason to pull a ~100 MB gmsh artifact into its environment — or into a documentation build.

The helpers

mfh_states(mat, dh, cv)one initial_state per quadrature point, as the nested [cell][qp] vector Ferrite loops expect
mfh_element!(Ke, re, cv, mat, ue, states, states_old, Δt; cache)element stiffness and internal force, small strain, plane strain
mfh_poro_element!(Ke, re, cvu, cvp, mat, ue, pe, …, Δt, mobility; u_range, p_range, cache)the coupled element — the equations
annulus_grid(Ri, Ro, nr, nθ)a structured annular sector, built by bending a rectangle — no gmsh
cylinder_sector_grid(Ri, Ro, H, nr, nθ, nz; grading)its 3-D twin, with geometric radial layers for a well problem

Keep two state arrays and swap only once a step has converged; material_response never mutates its argument, so a rejected Newton iteration is undone by not swapping.

julia
states     = EXT.mfh_states(mat, dh, cv)
states_old = EXT.mfh_states(mat, dh, cv)
cache      = MaterialCache()

for cell in CellIterator(dh)
    reinit!(cv, cell)
    fill!(Ke, 0); fill!(re, 0)
    ed = celldofs(cell)
    EXT.mfh_element!(Ke, re, cv, mat, u[ed],
                     states[cellid(cell)], states_old[cellid(cell)], Δt; cache)
    assemble!(assembler, ed, Ke, re)
end

A worked mechanical model is the thick-walled cylinder; the coupled element drives the ARMA 2011 well test.

The two-field element takes the mobility as an argument, one value per quadrature point: the permeability follows the apertures through a self-consistent solve, so the driver evaluates it once per step from the converged state and the scheme stays implicit in and explicit in .

julia
for c in eachindex(states), q in eachindex(states[c])
    K = transport_property(mat, states[c][q])
    mob[c][q] = K === nothing ? zero(SymmetricTensor{2,3}) : to_tensors(K) / μ
end

Two traps

A DofHandler numbers dofs by cell traversal, not node order — read nodal results with evaluate_at_grid_nodes(dh, u, :u), never by indexing u with a node number. And a MaterialCache is an unlocked Dict: use one per thread, or none.

Other codes

Only Ferrite is wired today. The contract itself is backend-agnostic, so a Gridap, FEniCSx or Abaqus UMAT driver needs the same two pieces — a per-point state container and an element routine — plus, outside Julia, a way to carry ε and σ across the boundary (voigt_strain and voigt_stress exist for that). Those are on the roadmap, not in the package.