General API

create_assembler(material, dh, ch; ah=nothing, quadrature_order=nothing, thickness=1.0)
create_assembler(materials_by_cellset, dh, ch; ah=nothing, quadrature_order=nothing, thickness=1.0)

Build a GenericMaterialAssembler from either a single material applied to all cells, or a dictionary mapping Ferrite cellset names to materials. Create the assembler once, after the DofHandler and ConstraintHandler are closed, and reuse it for the whole simulation.

Arguments

• material / materials_by_cellset: a single AbstractMaterial or a Dict{String,<:AbstractMaterial}. The dictionary form assigns different materials to different cellsets of the same grid, for example a stiff inclusion embedded in a soft matrix. The named cellsets must cover the active cells of dh exactly once. Cells outside get_cells(dh) are ignored.
• dh::DofHandler: a closed Ferrite DofHandler. Every sub-dof-handler of a mixed-element grid is supported.
• ch::ConstraintHandler: a closed Ferrite ConstraintHandler. Its Dirichlet pattern allocates both the linear and tangent sparse matrices.
• ah=nothing: optional user-supplied object. When it is non-nothing, the assembler calls create_alpha_values for a per-CellValues alpha-values object, then calls reinit!(av, cellid) for each cell. Material element routines apply alpha_value(av, qp) in the integration weight of stiffness and residual contributions, and compute_stresses applies the same factor to stress output. With ah === nothing, the default alpha_value(::Nothing, ::Int) = 1.0 applies.
• quadrature_order=nothing: optional quadrature order. An Int sets one order for every sub-dof-handler, a Vector{Int} sets one order per sub-dof-handler in order, and nothing auto-detects from the field interpolation (typically order(ip) + 1 for order > 1, else 2). The order is keyed by sub-dof-handler, not by material, so two materials sharing one SubDofHandler also share its quadrature rule. To give the cellsets of a multi-material grid different orders, declare one SubDofHandler per cellset and pass a Vector{Int}.
• thickness=1.0: scalar multiplier applied to the element stiffness and residual as ke .*= thickness; re .*= thickness. Omit it for 3D problems, where the default 1.0 is a no-op. For 2D plane problems, pass the physical out-of-plane thickness the model requires, for example thickness=0.1 for a millimetre-thick sheet. The wrappers PlaneStrain and PlaneStress require no thickness; their stress is per-unit-thickness by convention.

Additional call shapes

Four positional call shapes are accepted, taking either a single material or a material dictionary, with or without thickness:

create_assembler(material, dh, ah, ch, quadrature_order)
create_assembler(material, dh, ah, ch, quadrature_order, thickness)
create_assembler(material_dict, dh, ah, ch, quadrature_order)
create_assembler(material_dict, dh, ah, ch, quadrature_order, thickness)

These forward to the keyword form.

What is precomputed

• Sparse-matrix allocation matching the Dirichlet pattern of ch
• Per-cell mapping from cellid to material and to sub-dof-handler
• One CellValues per sub-dof-handler
• The linear stiffness of materials with is_linear == true (Hooke, Hooke2D), preassembled once into K_linear
• One AbstractMaterialState per owned nonlinear cell per quadrature point
• Per-owned-cell result buffers used by the two-phase assembly
• A stress cache of shape (max_nquadpoints, ncells)

See GenericMaterialAssembler for the full list of fields.

GenericMaterialAssembler

Stateful, preallocated Ferrite-based assembler for material models implementing the FerriteSolidMechanics material model interface.

The assembler owns the preassembled linear stiffness, the current nonlinear tangent and residual buffers, the material mapping for each active cell, and the quadrature point state vectors for nonlinear materials. It also caches sub-dofhandler metadata, task-owned workspaces for parallel assembly, per-cell element results used during scatter, and reusable stress-output buffers.

Construct it with create_assembler, then call stiffness_matrix, compute_forces, compute_stresses, and update_states!. The Concepts page describes the assembly and state-management model.

stiffness_matrix(assembler, u; dt=0.0)

Assemble and return (K, r), the current tangent stiffness matrix and the internal residual vector for the displacement vector u.

K is the algorithmic tangent. For materials that obtain their tangent via automatic differentiation (e.g. NeoHooke, VEPD_Detrez2010, VEVP_Zhao2021_AD, VEVP_MOAMMM) it is exact up to AD numerics, and for materials with an analytic tangent it is the linearized form returned by the model.

dt is the time step size of the current load step and is forwarded to every material's element routine. Rate-independent models (Hooke, Hooke2D, NeoHooke, J2Plasticity) typically ignore it, while viscoelastic and viscoplastic models use it to evolve internal variables. For viscous models, dt=0.0 freezes the viscous evolution rather than producing rate-independent behaviour.

Assembly runs in two passes:

1. Each owned nonlinear cell computes its element ke and re in parallel (OhMyThreads.@tasks, greedy scheduler)
2. Per-cell results are scattered into the global matrix and vector, and reduced across MPI ranks when MPI.COMM_WORLD has more than one process

Linear cells, those whose material has is_linear == true, are preassembled once when the assembler is created and are added to the tangent on every call without re-evaluation.

Under MPI this is a collective call that every rank must reach. A local failure on one rank is made collective before it propagates, so peer ranks throw instead of blocking in the reduction. Use try_stiffness_matrix when the outer solver should recover from local constitutive failures rather than abort.

K and r are the assembler's reusable buffers and are overwritten by the next assembly call. copy them to keep a previous tangent or residual.

try_stiffness_matrix(assembler, u; dt=0.0)

Non-throwing wrapper around stiffness_matrix for adaptive outer solvers that reject the current load or time step and retry with a smaller dt when a local constitutive update fails. Returns a named tuple:

• (converged = true, K = K, r = r, error = nothing) on success
• (converged = false, K = nothing, r = nothing, error = err) when a recoverable LocalAssemblyFailure occurs

Only LocalAssemblyFailure is treated as recoverable. Every other exception is rethrown, including a MethodError in a material implementation and a generic DomainError that does not subtype LocalAssemblyFailure. A material whose local constitutive failure should be retried by an adaptive outer solver must therefore throw an exception subtyping LocalAssemblyFailure. Threaded assembly may wrap failures in task exceptions; a wrapped failure is recoverable only when every wrapped leaf exception is a LocalAssemblyFailure.

On failure, reject the current outer step, restore the quadrature states with revert_states!, restore the displacement vector to the last accepted step, reduce dt, and assemble again. Do not call update_states! for a failed attempt; commit only after the outer Newton solve has converged. The assembler's matrix and residual buffers are meaningless after a failed attempt, so use only the returned error.

In MPI runs, unrecoverable local errors are made collective before they are thrown, so peer ranks do not wait in a later reduction. Recoverable failure flags are also reduced across MPI.COMM_WORLD before the global tangent and residual reductions. If one rank fails recoverably, every rank receives converged = false, and ranks without the original exception receive RemoteAssemblyFailure as error.

compute_forces(assembler, u; dt=0.0)

Return only the internal force (residual) vector for displacement u.

compute_forces wraps stiffness_matrix and discards the tangent matrix. It runs the same assembly path, so the per-cell element loops still compute both ke and re. Inside a Newton solve, call stiffness_matrix directly and ignore the matrix instead, which avoids a second assembly.

dt has the same meaning as in stiffness_matrix.

The returned vector is the assembler's reusable residual buffer and is overwritten by the next assembly call. copy it to keep a previous residual.

compute_stresses(assembler, u; dt=0.0)

Evaluate quadrature point stresses for the current displacement vector u and return them as an array of Tensor values, specifically Tensor{2,dim}.

The returned array has shape (max_nquadpoints, getncells(grid)) and is indexed σ[qp, cellid] by Ferrite cell id. It is sized for the whole grid, so when the DofHandler covers only a subset of cells, the columns of inactive cells stay zero and the active cells keep their own ids as column indices; they are not compacted to 1:n_active. For mixed sub-dof-handlers, entries above a cell's own quadrature count remain zero. The array is the assembler's reusable internal stress cache, and a later compute_stresses call on the same assembler overwrites it. copy it to keep a previous stress field.

Each quadrature point stress is multiplied by alpha_value(av, qp) when the assembler was built with a non-nothing ah keyword, and by 1.0 otherwise. See alpha_value.

dt defaults to 0.0 and is passed to the material's stress routine for API consistency. For postprocessing, call compute_stresses(assembler, u) and leave that default in place. Stress postprocessing is expected to be non-evolving: bundled stateful materials report the current trial or the committed stress without advancing history or committing state. compute_stresses is not a substitute for stiffness_matrix, compute_forces, or update_states!. Perform the time step update with assembly, then call update_states! after the outer solve has converged.

compute_stresses loops over all owned cells, linear and nonlinear, and runs in parallel with OhMyThreads.@tasks. When MPI.COMM_WORLD has more than one process the per-rank results are summed across ranks, so every rank ends up with the same global stress field. It is therefore a collective call that every rank must reach. A local failure on one rank is made collective before it propagates, so peer ranks throw instead of blocking in the reduction.

update_states!(assembler)

Commit the trial states of all owned nonlinear quadrature points after a converged Newton step.

update_states! calls update_state! on every quadrature point state of every owned nonlinear cell, copying the trial values into the committed slot. Call it once per load step, after Newton has converged, and never inside the Newton loop. If the step is rejected instead, call revert_states! to discard the trial values.

update_state! either mutates the state in place and returns nothing, or returns a replacement AbstractMaterialState. A returned state is stored back into the quadrature point state vector.

Cells with a linear material hold no history and are skipped, as are cells owned by other MPI ranks.

revert_states!(assembler)

Discard the trial states of all owned nonlinear quadrature points and restore their last committed values.

revert_states! calls revert_state! on every quadrature point state of every owned nonlinear cell, copying the committed values back into the trial slot. Call it when a trial step is rejected, typically before retrying with a smaller time step.

revert_state! either mutates the state in place and returns nothing, or returns a replacement AbstractMaterialState. A returned state is stored back into the quadrature point state vector.

Cells with a linear material hold no history and are skipped, as are cells owned by other MPI ranks.

recommended_solve_settings(dh::DofHandler; cores_per_node=Sys.CPU_THREADS,
                           memory_per_node_gb=nothing, mumps_available=false,
                           assembly_bound=false) -> NamedTuple

Return the launch settings for a run on one node: which linear solver to use, how many MPI ranks to start, how many Julia threads to give each, and how many BLAS threads.

There are two ways to solve K \ r each Newton iteration, and the choice is worth up to 5× here:

• K \ r (:replicated_lu) – every rank solves the whole system by itself, so each pays the full memory of one solve
• distributed_solve (:distributed_solve) – the ranks split one solve between them, so each pays a share of the memory

distributed_solve needs import MUMPS; pass mumps_available=true once you have it.

recommended_solve_settings returns (; solver, ranks_per_node, threads_per_rank, blas_threads, memory_gb, fits, note). memory_gb is what one rank needs for K \ r, from estimated_replicated_memory_gb. fits is true when all ranks together stay within 85% of memory_per_node_gb, budgeting 2.6× memory_gb per rank. That budget is deliberately pessimistic and reports false before a run would actually fail; Memory explains the headroom behind the 2.6×. VEVP_Zhao2021_AT measured 9× the estimate and exceeds even this budget, so size a model with many Newton iterations per load step from its own measurement. fits is missing if memory_per_node_gb is not specified or if distributed_solve has been selected: In this case, the factorization is split across the ranks, so the replicated-copies budget checked by fits does not apply. In 2D above a few million dofs also check that all ranks' copies of K fit in the node's memory, no matter the solver: every rank holds its own full copy, and at 16M dofs 18 such copies alone exceed the memory of a 256 GB node.

Pass assembly_bound=true for a 2D run that spends more time building element matrices than solving, which is typical of rate-dependent models that run their own local Newton solve per quadrature point. The flag moves the 2D layout toward more ranks. Raise the rank count further yourself if assembly dominates and the memory allows it; Rank and thread layout shows where the highest rank count stops fitting as the mesh grows.

s = recommended_solve_settings(dh; cores_per_node=72, memory_per_node_gb=256,
                               mumps_available=true)
LinearAlgebra.BLAS.set_num_threads(s.blas_threads)
# launch with s.ranks_per_node ranks and s.threads_per_rank threads each

The rules below were measured on the 72-core nodes of NHR@FAU Fritz:

recommended_solve_settings holds the default 2D rank count at cores ÷ 8 and, with mumps_available=true and memory_per_node_gb given, returns distributed_solve once that many copies would exceed 85% of the node's memory. The rank and thread counts come from one mesh per case, so treat them as starting points and better time your own problem before committing to a large run. Higher-order elements were not measured; at equal dof count they change the sparsity of K, so remeasure the thresholds and the memory behind fits before relying on them for another element order. Choosing the linear solver states the timings behind these rules, and Reproducing these on your cluster shows how to check them on your own cluster.

recommended_blas_threads(dh::DofHandler; ncores=Sys.CPU_THREADS, ranks_per_node=nothing) -> Int

Return the number of BLAS threads for the K \ r solve the driver runs each Newton iteration. Set it in the driver before solving: LinearAlgebra.BLAS.set_num_threads(recommended_blas_threads(dh)).

recommended_blas_threads returns ncores for a 3D grid whose rank has the node to itself, and 1 in 2D or when several ranks share a node. ranks_per_node defaults to the number of ranks on this rank's node, read from MPI. Julia starts BLAS at Sys.CPU_THREADS ÷ 2, so the returned count is higher in the 3D single-rank case and lower in the other two.

BLAS threading pays off in 3D when the rank owns the node: at 108k dofs the solve took 312 s at 1 BLAS thread against 145 s at 36. In 2D it measured slower at every thread count. Ranks sharing a node each start their own BLAS threads, which cost distributed_solve 6.3× per solve at 18 ranks.

recommended_solve_settings returns all launch settings at once. BLAS threads states the measurements behind this rule, including the mesh size below which BLAS threading stops paying.

estimated_replicated_memory_gb(ndofs, dim) -> Float64

Estimate the memory one MPI rank needs, in GB, when the driver solves K \ r directly.

Solving K \ r builds a decomposition of K that is much bigger than K itself, and every rank builds its own copy. That copy is what fills up a node: running 18 ranks needs 18 times this number. distributed_solve splits one copy across the ranks instead, so it needs far less memory per rank.

Memory grows faster than the problem: ndofs^1.43 in 3D and ndofs^1.09 in 2D. Doubling a 3D mesh's dofs therefore costs about 2.7× the memory, so a size that fits comfortably can stop fitting after one refinement.

The fitted curve is scaled to NeoHooke over three load steps, on linear hexahedra in 3D and linear quadrilaterals in 2D, and the predicted peak memory stays within 15% of every point measured on the 72-core nodes of NHR@FAU Fritz. Higher-order elements were not measured; at equal dof count they change the sparsity of K, so remeasure before relying on the estimate for another element order. An expensive material and a longer run both raise the peak above it, together by 2.6×, which is the headroom recommended_solve_settings budgets in its fits field. Either the expensive material or the longer run reaches that budget on its own, and VEVP_Zhao2021_AT exceeded it at 9× the estimate, so take the memory of a long run, or of a model with many Newton iterations per load step, from its own measurement.

Memory documents the range of grid sizes the fit was validated over, the measured series behind the 2.6×, and the size above which each rank's own copy of K rather than the factorization sets the limit.

distributed_solve(K, residual, comm; verify_rtol=1e-8, workspace_pct=nothing) -> du

Solve K * du = residual with the distributed sparse direct solver MUMPS, returning the update du on every rank. It is a drop-in replacement for du = K \ residual in an MPI Newton loop:

u .-= distributed_solve(K, residual, MPI.COMM_WORLD)

The method comes from a package extension. Add MUMPS.jl to the environment and load it with import MUMPS, which is what defines the method. Calling distributed_solve beforehand raises a MethodError.

Every rank passes its own copy of the assembled K and residual. MUMPS factorizes across comm and broadcasts the solution back, so all ranks return the same du. residual is left unchanged.

Before returning, distributed_solve checks du against K * du ≈ residual to a relative tolerance of verify_rtol and throws a DistributedSolveError when that check fails. The check exists because MUMPS reported success on a cluster run whose solution was wrong. It costs one sparse matrix-vector product per call (one per Newton iteration; on each rank's own copy and without synchronizing), which is negligible next to the factorization. A working solve satisfies it to roughly machine precision. Set verify_rtol=0 to skip it.

MUMPS sizes its internal workspace from an estimate, and fails the factorization when that estimate turns out too small (INFO(1) = -9), the most frequent failure on the problems measured here. distributed_solve retries that case automatically with progressively more workspace, so that a solve claims extra memory and completes instead of ending the run. workspace_pct sets that percentage (MUMPS ICNTL(14)) explicitly and turns the retries off, which caps the memory a single solve may claim.

Use this in 3D above between roughly 4k..20k dofs. With K \ residual every rank builds the whole factorization, so time and memory both grow with the rank count, while MUMPS splits one factorization across the ranks instead. Measured on the 72-core nodes of NHR@FAU Fritz against the best replicated configuration (one rank × 72 cores, BLAS 36), on 3D hexahedral NeoHooke problems, where both methods produce the same solution to within 1e-9 relative:

The two ways of solving cost the same per solve at about 3k dofs, and the 2.2k row covers eight repeats each whose ranges overlap, so the replicated solve is the more straightforward choice in the vicinity of the crossover. In 2D, K \ residual was faster at every size measured (291k: 52 s vs 62 s; 4.0M with Hooke: 30 s vs 39 s), both at BLAS 1, so use it there with fewer ranks than the 18 the distributed layout runs, since each rank's whole factorization copy must fit in node RAM. At 16M 2D dofs one replicated rank completed using 95 GB of memory while 18 ranks running distributed_solve were killed on a 256 GB node: each of them holds its own 4.3 GB copy of K on top of its share of the factorization, so at that scale fewer ranks save more memory than splitting the factorization does. recommended_solve_settings applies these rules to a given grid, keeping 18 ranks × 4 threads in 3D rather than the 36 × 2 layout the 273k row used.

Launch it the opposite way from the threaded K \ residual, which runs one rank with all 72 cores as threads: distributed_solve takes its parallelism from MPI processes instead, so run many ranks per node (18 to 36 on 72 cores) with a single BLAS thread each. At 18 ranks, 4 BLAS threads per rank measured 6.3× slower per solve than 1 thread. recommended_blas_threads applies to the replicated solve. The factorization communicates heavily between ranks, so these measurements all used the cluster's own MPI module rather than the build MPI.jl installs by default. Use one node while the problem fits its RAM. At 273k dofs the run took 102 s on 1 node against 123 s on 2 and 76 s on 4, all at 18 ranks per node, while 36 ranks on the single node took 62 s and beat every one of them. At 985k dofs it took 377 s on 2 nodes, 708 s on 4 and 821 s on 8, the factorization itself accounting for 122 s, 233 s and 270 s per solve, so adding nodes bought no solve time at either size. Element assembly does scale with the ranks, from 0.168 s to 0.064 s per Newton iteration at 273k dofs over 18 to 72 ranks, but that saving is negligible next to the factorization. Furthremore, more nodes help only once the memory requires them. The MPI tutorial covers layout, mesh-size, and memory scaling.

DistributedSolveError(phase, infog1, [relative_residual]) <: LocalAssemblyFailure

distributed_solve did not return a usable du.

phase is :factorize or :solve when MUMPS reported the failure itself, and infog1 is then its INFOG(1)/INFO(1) error code (-10 singular, -13 allocation failure, -9 workarray too small after every workspace retry has been exhausted). phase is :verify when MUMPS reported success but K * du did not reproduce residual, and relative_residual is then how far apart the two were, relative to residual, while infog1 is 0.

DistributedSolveError subtypes LocalAssemblyFailure, so a driver that shrinks dt and retries the step can catch it.

LoadHandler(assembler; quadrature_order=nothing)
LoadHandler(dh; thickness=1.0, quadrature_order=nothing)

External load container, built and used the way Ferrite's ConstraintHandler is:

ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> [0.0, 0.0], [1, 2]))
close!(ch)

lh = LoadHandler(assembler)
add!(lh, Traction("top", (x, t) -> Vec(0.0, -p * t)))
add!(lh, NodalForce("tip", x -> Vec(0.0, -F); distribute=true))
close!(lh)

external_forces! then returns the assembled external force vector for a given t. The handler follows the pattern and keeps the same name as the LoadHandler of FerriteAssembly.jl (K. A. Meyer).

Constructing from a GenericMaterialAssembler takes the DofHandler, the thickness and the quadrature orders from it, so external and internal forces are integrated consistently. Passing a DofHandler instead is for drivers that assemble no internal forces; give the same thickness the assembler uses, because facet and body integrals scale with it exactly as ke and re do.

quadrature_order accepts an Int for one order across all sub-dof-handlers, a Vector{Int} for one order each, and nothing to take the assembler's orders or, without an assembler, the interpolation default.

Adding loads after closing

close! prepares the entries added since the last call and leaves earlier ones untouched, so loads may be added part-way through a simulation and the handler closed again:

add!(lh, Traction("second_stage", (x, t) -> Vec(0.0, -q * (t - t_switch))))
close!(lh)

Changing the magnitude, direction or profile of an existing load needs no such call, because the load function is evaluated on every external_forces! call and may close over computed results:

stage = Ref(1)
add!(lh, Traction("top", (x, t) -> stage[] == 1 ? Vec(0.0, -p) : zero(Vec{2})))

The same pattern removes a load, since a function returning zero contributes nothing to the external force vector.

Replication under MPI

LoadHandler assembles serially and, under MPI, redundantly: every rank builds the complete external force vector and no reduction is performed. This matches the replicated-data design of the assembler, whose residual is already complete on every rank after its own Allreduce. Every rank must add the same loads in the same order so that the per-rank copies agree. An MPI.Allreduce over the external force vector would sum those identical copies and apply the load nranks times.

external_forces!(lh::LoadHandler, t=0.0)

Assemble the external force vector of every load in lh at t and return it.

t is passed unchanged as the second argument of each load function. The driver chooses its meaning: physical time for rate-dependent materials, or a load factor for quasi-static ramping. Pass the same value to Ferrite's update!(ch, t) so that prescribed displacements and loads refer to the same instant.

Dead loads do not depend on the displacement, so call external_forces! once per time step, outside the Newton loop, and subtract the result from the internal force returned by stiffness_matrix:

f_ext = external_forces!(lh, t)
for iter in 1:max_newton
    K, f_int = stiffness_matrix(assembler, u; dt=dt)
    residual .= f_int .- f_ext
    apply_zero!(K, residual, ch)
    norm(residual) <= tol * max(norm(f_ext), 1.0) && break
    u .-= K \ residual
end

External loads contribute nothing to K.

Under MPI every rank assembles the complete vector and no reduction is performed; see LoadHandler.

The returned vector belongs to lh. Calling external_forces! again overwrites it, so copy the result when a step's values are still needed afterwards.

Throws when loads have been added since the last close!(lh), which must be called before the added loads can be assembled.

Traction(set, f)

Surface traction on a facet set, in force per unit area in 3D and force per unit length in 2D.

set is a facetset name or a set of FacetIndex. A facet of a 2D cell is an edge, so the same type covers area loads in 3D and edge loads in 2D.

f returns the traction vector and may be written in any of three forms, resolved once when the LoadHandler is closed:

Traction("top", x -> Vec(0.0, -p))                   # f(x)
Traction("top", (x, t) -> Vec(0.0, -p * t))          # f(x, t)
Traction("top", (x, t, n) -> -p * t * n)             # f(x, t, n)

x is the spatial coordinate of the quadrature point, t is the value passed to external_forces!, and n is the outward facet normal. Because f receives x, any spatial profile is expressed directly, for example (x, t) -> Vec(0.0, -p * (x[1] / L) * t) for a load growing linearly along the x-axis.

n is evaluated on the undeformed geometry, so a normal-aligned traction is a dead load, fixed in direction at the start and not rotating as the surface deforms. Follower loads, whose direction tracks the deforming surface, would contribute an extra term to the tangent and are not available. For the common normal-aligned case, Pressure carries the sign convention.

Pressure(set, f)

Normal pressure on a facet set, in the same units as Traction.

f returns a scalar and accepts the forms f(x), f(x, t) and f(x, t, n). The applied traction is -f * n with n the outward facet normal, so a positive value acts inward and compresses the body.

Pressure("hole", (x, t) -> p * t)

Pressure is equivalent to Traction(set, (x, t, n) -> -f(x, t) * n) and exists to keep that sign convention in one place. The normal is taken on the undeformed geometry, so this is a dead load; see Traction.

NodalForce(set, f; distribute=false)

Point force applied at the nodes of a node set, in force.

set is a nodeset name or a collection of node ids. f returns the force vector and accepts the forms f(x) and f(x, t), where x is the node coordinate. The facet normal is not defined at a node, so the f(x, t, n) form is not available here.

distribute=false, the default, applies the returned value at every node of the set. distribute=true treats the returned value as the resultant of the whole set and splits it equally over the nodes, so the sum over the set is the value itself.

NodalForce("tip", x -> Vec(0.0, -F))                      # F at each node
NodalForce("tip", x -> Vec(0.0, -F); distribute=true)     # F in total

Equal splitting matches the finite element load distribution only for first-order shape functions. Above first order the shape function integration no longer splits the load equally, so corner and midside nodes take different shares. To spread a resultant over a surface on such a mesh, use Traction with magnitude F / A instead.

Only nodes carrying degrees of freedom can receive a force. Closing the handler throws when the set contains nodes without dofs.

BodyForce(f)
BodyForce(set, f)

Body force over cells, in force per unit volume in 3D and force per unit area in 2D.

The one-argument form applies to every cell carrying degrees of freedom. The two-argument form takes a cellset name or a collection of cell ids. f returns the force vector and accepts the forms f(x) and f(x, t), where x is the spatial coordinate of the quadrature point.

BodyForce(x -> Vec(0.0, -rho * g))                        # self-weight
BodyForce("rotor", (x, t) -> rho * omega(t)^2 * x)        # centrifugal, per cellset

AbstractLoad

Supertype of the external load specifications Traction, Pressure, NodalForce and BodyForce.

Loads are Neumann boundary conditions: they contribute to the external force vector and never to the tangent. Prescribed displacements are Dirichlet boundary conditions and are not part of this hierarchy; use Ferrite's Dirichlet and ConstraintHandler for those.

Each load type is named after the units of the value its function returns. The 2D and 3D units differ only because the entity integrated over does.

TimeStepController(; dt_min=NaN, dt_max=NaN, shrink=0.5, grow=2.0,
                   max_rejections=typemax(Int))

Adaptive time step controller for user-owned load and time loops.

The controller runs no assembly, Newton iterations, rollback, or MPI communication. It only updates scalar time step sizes after the caller has accepted or rejected a trial step, which makes it usable with any material model and with serial, threaded, or MPI assembly.

dt_min and dt_max are optional bounds, and NaN disables the corresponding lower or upper bound. shrink must satisfy 0 < shrink < 1. grow must be finite and satisfy grow >= 1. Both factors are fixed when the controller is constructed; accept_step! and reject_step! do not adapt them. max_rejections = N allows N consecutive rejected steps since the last accepted step and throws on rejection N + 1. The default leaves that count practically unbounded.

rejections counts the rejected steps since the last accepted one. reject_step! increments it, accept_step! resets it to zero, and reset_controller! clears it without returning a new step size. It is the only field meant to change after construction; leave dt_min, dt_max, shrink, grow and max_rejections at the values given to the constructor.

Use one TimeStepController per Newton loop. The mutating functions assume a single caller and are not safe to call from several tasks at once.

TimeStepControllerExhausted

Exception thrown by reject_step! when the controller has no smaller time step to offer, because the proposed step would violate dt_min or the consecutive rejection budget.

accept_step!(controller, dt)

Return the next trial time step after an accepted step.

accept_step! resets the controller's consecutive rejection counter, multiplies dt by controller.grow, and clamps the result to controller.dt_max when that upper bound is enabled. dt_max = NaN disables the upper clamp.

reject_step!(controller, dt)

Return the next trial time step after a rejected step.

reject_step! computes dt * controller.shrink, checks the lower bound and the rejection budget, then increments the controller's consecutive rejection counter on success. It throws when the new step would fall below controller.dt_min and that lower bound is enabled. It does not clamp rejected steps to dt_min. dt_min = NaN disables the lower-bound check. The controller state is left unchanged when reject_step! throws.

reset_controller!(controller)

Reset the controller's consecutive rejection counter to zero and return the controller.

AbstractMaterial

Supertype for all material models used with GenericMaterialAssembler.

AbstractMaterialState

Supertype for all quadrature point state objects that store internal history.

LocalAssemblyFailure

Abstract exception type for recoverable local assembly failures. Public non-throwing helpers such as try_stiffness_matrix catch this type and return it to the caller, so an outer load step controller can reduce the step size and retry.

Those helpers catch only exceptions that opt in by subtyping LocalAssemblyFailure. Exceptions from invalid inputs or buggy material methods are rethrown.

NoState

Placeholder state for materials that do not require internal history.

material_response(material, ε_or_F, state, dt, cache=nothing) -> (stress, tangent, new_state)

Constitutive update at a single quadrature point. A typical material implements material_response, kinematics, and, for a history-dependent model, a state type.

Arguments

• ε_or_F: strain measure matching kinematics(material) – ε::SymmetricTensor{2,dim} for SmallStrain, F::Tensor{2,dim} for FiniteStrain
• state: the quadrature point state object, NoState() for a stateless material. Read the committed history from it, by convention the previous-side fields. Do not commit.
• dt: time increment of the current step, where 0.0 freezes viscous flow
• cache: scratch object from allocate_material_cache, or nothing

Returns

• stress: σ for SmallStrain, P for FiniteStrain
• tangent: ∂σ/∂ε for SmallStrain, ∂P/∂F for FiniteStrain
• new_state: the trial state produced by this evaluation

The generic element routine accepts three forms of new_state:

1. the passed state itself, for a stateless material or a material that wrote its trial fields in place
2. a snapshot of internal variables that is not an AbstractMaterialState, stored via set_trial!(state, new_state)
3. a new AbstractMaterialState, which replaces the quadrature point state slot and must have the same concrete type as the state it replaces

The generic element routine pins one state type per cell, taken from the first quadrature point.

material_response must be AD-transparent where the material relies on automatic differentiation. It must throw a LocalAssemblyFailure subtype for a recoverable local-solve failure.

The material_response name and the AbstractMaterial/AbstractMaterialState pair come from MaterialModels.jl by K. Auth and contributors. K. A. Meyer's MaterialModelsBase.jl was developed from that package as an implementation-independent interface. The design here follows the MaterialModelsBase conventions: a single response function returning (stress, tangent, new_state), strain-measure-driven stress/tangent pairs, and the cache hook. Materials implementing the MaterialModelsBase interface can be used directly via FromMaterialModelsBase.

material_stress(material, ε_or_F, state, dt, cache=nothing) -> stress

Return the postprocessing stress at one quadrature point, without the tangent.

material_stress is the optional stress-only companion to material_response. The generic _compute_stress_qp fallback calls it to obtain the output stress: Cauchy σ for FiniteStrain, small-strain σ for SmallStrain. The default method computes that stress from material_response at the passed state and discards the tangent and the trial state, so a plain constitutive material needs no method of its own.

Define a method when computing the tangent is expensive, or when the output stress follows a model-specific convention. An AD-based history model, for example, can return the stress directly from its local solve and skip the automatic differentiation. Override material_stress in preference to the lower-level _compute_stress_qp: material_stress receives the strain measure ε or F directly, whereas _compute_stress_qp receives the cellvalues and the quadrature point index and leaves the strain extraction to the method.

A material_stress method must be deterministic and must not commit or evolve state. A rate-dependent model, for example, passes dt = 0.0 to its local solve by default, so it reports the stress at the current time and state without evolving state.

kinematics(material) -> AbstractKinematics

Declare the strain measure of a material: SmallStrain or FiniteStrain.

Extension point. Every material that uses the generic element routine, meaning every model that implements material_response, must define kinematics. Only a material that provides solely its own _assemble_element! method, and is not meant to be wrapped, may omit it.

AbstractKinematics

Trait supertype describing which strain measure a material consumes. kinematics returns it, and the generic element routine uses it to select what to evaluate at each quadrature point and which weak form to integrate. Concrete traits: SmallStrain, FiniteStrain.

SmallStrain <: AbstractKinematics

Material works on the small-strain tensor ε = sym(∇u), a SymmetricTensor{2,dim}, and returns the Cauchy stress σ and the algorithmic tangent D = ∂σ/∂ε.

FiniteStrain <: AbstractKinematics

Material works on the deformation gradient F = I + ∇u, a Tensor{2,dim}, and returns the first Piola–Kirchhoff stress P and the tangent ∂P/∂F.

tangent_symmetry(material) -> AbstractTangentSymmetry

Declare whether a SmallStrain material's tangent D = ∂σ/∂ε is major-symmetric. Defaults to Unsymmetric.

Returning MajorSymmetric halves the tangent-integration work, worth roughly 1.4× on small-strain element assembly. The declaration is a promise about the values, not about the type. SymmetricTensor{4} carries only the minor symmetries D[i,j,k,l] == D[j,i,k,l] == D[i,j,l,k], which do not imply major symmetry. Check before claiming it:

D ≈ permutedims(D, (3, 4, 1, 2))

FiniteStrain materials do not use this trait, because their ∂P/∂F is never assumed symmetric.

AbstractTangentSymmetry

Trait supertype describing whether a SmallStrain material's tangent D = ∂σ/∂ε is major-symmetric. tangent_symmetry returns it, and the generic element routine uses it to select the tangent-integration loop. Concrete traits: MajorSymmetric, Unsymmetric.

MajorSymmetric <: AbstractTangentSymmetry

D[i,j,k,l] == D[k,l,i,j], so the element stiffness is symmetric. The element routine integrates only its lower triangle and mirrors the upper triangle.

Declaring this trait for a tangent that is not major-symmetric silently discards the antisymmetric part and yields a wrong element stiffness, so Newton converges slowly or not at all. The residual stays correct, which makes the error easy to miss.

Unsymmetric <: AbstractTangentSymmetry

No symmetry is assumed, so the element routine integrates every entry of ke. This is the default, and the correct choice for non-associated flow, nonlinear kinematic hardening, and anisotropic damage.

set_trial!(state, new_state)

Store the trial result new_state returned by material_response into the quadrature point state without committing it. The default method covers the current/previous layout:

set_trial!(state, new) = (state.current = new; nothing)

Override it for a state with a different layout; VEPD_Detrez2010 is a field-wise example. The generic element routine does not call set_trial! when material_response returns the state itself or a replacement AbstractMaterialState.

allocate_material_cache(material) -> cache

Allocate a reusable scratch object that is passed as the cache argument of material_response. The default method returns nothing. Element assembly calls allocate_material_cache once per element, and the generic stress fallback calls it once per quadrature point.

Override it to reuse buffers across the quadrature points of an element, such as workspace arrays for a local Newton solve. A material whose cache is expensive to build should also override material_stress.

is_linear(material)

Return true when material has a displacement-independent tangent stiffness. GenericMaterialAssembler preassembles linear materials once.

create_state(material)

Create the quadrature point state object of material. A stateless material uses the default method, which returns NoState().

update_state!(state)

Commit one converged quadrature point state at the end of a load or time step.

update_state! either mutates a mutable state in place and returns nothing, or returns a new AbstractMaterialState for an immutable state. update_states! stores a returned state back into the assembler state vector.

revert_state!(state)

Restore one quadrature point state to the previously committed value.

revert_state! either mutates a mutable state in place and returns nothing, or returns a new state value for an immutable state. A caller that owns the state container is responsible for storing that returned value.

copy_state!(dest, src)

Copy fields from src to dest recursively, using reflection. copy_state! copies arrays in place where possible and returns immutable values, so a caller that owns the storage can replace the destination value. This is the generic history-state fallback for arbitrary structs.

AbstractHyperelastic <: AbstractMaterial

Supertype for stateless hyperelastic materials defined by a strain energy density Ψ(C). A subtype implements the energy function and receives the complete material interface:

struct MyModel <: AbstractHyperelastic
    μ::Float64
end
FerriteSolidMechanics.Ψ(C, mp::MyModel) = 0.5 * mp.μ * (tr(C) - 3)  # + volumetric part

Derived defaults, each overridable with an analytic expression:

• constitutive_driver(C, mp) – S = 2∂Ψ/∂C and ∂S/∂C = 2∂²Ψ/∂C² via Tensors.hessian
• _second_piola(C, mp) – S = 2∂Ψ/∂C via Tensors.gradient
• material_response, kinematics (= FiniteStrain()), _compute_stress_qp, compute_PK1_3D and update_state_from_3D!, so PlaneStrain and PlaneStress wrapping works without additional methods

Ψ(C, material)

Strain energy density as a function of the right Cauchy–Green tensor C = FᵀF. This is the extension point for AbstractHyperelastic materials. Extend it as FerriteSolidMechanics.Ψ.

FromMaterialModelsBase(material; kinematics=SmallStrain())

Wrap a material implementing the MaterialModelsBase.jl interface, for example a model from MechanicalMaterialModels.jl, for use everywhere a FerriteSolidMechanics AbstractMaterial is accepted: create_assembler, stiffness_matrix, compute_stresses, update_states!, try_stiffness_matrix, and so on.

The kinematics keyword declares which strain measure the wrapped material consumes, since MaterialModelsBase selects it by dispatch on the input type:

• SmallStrain, the default, passes ε::SymmetricTensor{2,dim} and expects (σ, ∂σ∂ε, new_state)
• FiniteStrain passes F::Tensor{2,dim} and expects (P, ∂P∂F, new_state)

The constructor is always available, while a package extension provides the interface methods behind it. Load the weak dependency first:

import MaterialModelsBase                    # activates the extension
using MechanicalMaterialModels: Plastic, LinearElastic, Voce,
                                ArmstrongFrederick, NortonOverstress

mat = FromMaterialModelsBase(Plastic(;
    elastic   = LinearElastic(E=210.0e3, ν=0.3),
    yield     = 100.0,
    isotropic = Voce(Hiso=50.0e3, κ∞=100.0),
    kinematic = ArmstrongFrederick(Hkin=200.0e3, β∞=200.0),
    overstress = NortonOverstress(tstar=1.0, nexp=2.0),
))
fem = create_assembler(mat, dh, ch)

Prefer import MaterialModelsBase over using MaterialModelsBase: both packages export a material_response function with the same semantics by design, so using both makes the unqualified name ambiguous.

Neither MaterialModelsBase nor MechanicalMaterialModels is registered in the General registry at the time of writing, so install them by URL: Pkg.add(url="https://github.com/KnutAM/MaterialModelsBase.jl"). MechanicalMaterialModels additionally depends on Newton.jl, which must be added first via Pkg.add(url="https://github.com/KnutAM/Newton.jl").

FromMaterialModelsBase stores the wrapped material's immutable state snapshots in an MMBState trial/commit pair, so update_states! and revert_states! work as for any bundled material. A local convergence failure, MaterialModelsBase.MaterialConvergenceError, is rethrown as MaterialModelsBaseConvergenceError, a LocalAssemblyFailure, so try_stiffness_matrix and TimeStepController-driven adaptive stepping work unchanged.

MMBState(current, previous)

Trial/commit state pair for FromMaterialModelsBase. Both fields hold immutable MaterialModelsBase state snapshots, and the default set_trial! semantics apply. update_state! commits previous ← current, and revert_state! rolls back current ← previous.

MaterialModelsBaseConvergenceError(error)

Recoverable LocalAssemblyFailure wrapping a MaterialModelsBase.MaterialConvergenceError thrown by a wrapped material's local solve, so that try_stiffness_matrix reports converged = false instead of aborting the outer solve.

_assemble_element!(ke, re, states, material, cellvalues, alphavalues, u, dt)

Element-level extension hook for assembling an element stiffness matrix ke and internal force vector re.

Most materials do not implement _assemble_element!. A generic fallback loops over the quadrature points, evaluates material_response, and integrates the weak form selected by kinematics. Implement _assemble_element! directly only when the element structure itself is material-specific, such as a hand-derived tangent assembled from nodal shape-gradient blocks (VEVP_Zhao2021_AT) or nonlocal regularization.

This hook does not reach mixed or multi-field formulations. The assembler builds one CellValues per SubDofHandler from its single field interpolation, so it rejects a SubDofHandler carrying more than one field when the assembler is created.

Extension point. The leading underscore indicates that users normally call higher-level assembler functions, but this method is exported so material packages can extend it. The Developer guide walks through a complete worked example.

A custom method is responsible for:

• looping over the quadrature points
• evaluating the strain or the deformation gradient
• computing the stress and the tangent
• multiplying by the integration weight, including the alpha_value(alphavalues, qp) * getdetJdV(cellvalues, qp) factor, which the assembler does not apply
• adding the contributions to ke and re
• writing trial values into the material-specific current/trial state fields, without committing them

The assembler calls update_state! only after the outer Newton step has converged. A custom method mutates ke, re, and the state in place and returns nothing; its return value is ignored.

_compute_stress_qp(material, cellvalues, alphavalues, qp, u_local, state, dt)

Material extension hook for postprocessing the stress at one quadrature point.

Extension point. The leading underscore indicates that users normally call higher-level assembler functions, but this method is exported so material packages can extend it.

A generic fallback derives the stress from material_stress, which itself defaults to material_response with the tangent and the trial state discarded. A material implementing the constitutive interface therefore needs no _compute_stress_qp method. For a tangent-free stress path, override the higher-level material_stress instead, which receives the strain measure directly. Define a _compute_stress_qp method only when the output convention needs the cellvalues and quadrature point context itself, as in element-structured models.

A _compute_stress_qp method returns the Cauchy stress at the quadrature point for a finite-strain material, or the Cauchy-like stress for a small-strain material. compute_stresses multiplies that value by alpha_value(alphavalues, qp) before storing it in the output array, so the material must not apply the alpha factor itself. A linear cell that owns no state is redispatched from state === nothing to NoState(). A material with neither a _compute_stress_qp method nor the material_response/kinematics pair throws an ArgumentError instead of returning a silent zero stress.

tangent_source(material) -> Function

Select where PlaneStrain and PlaneStress take the 3D stress and tangent of a FiniteStrain material from.

Defaults to material_response, whose tangent the wrappers restrict (plane strain) or condense (plane stress) directly. Returning compute_PK1_3D instead makes them differentiate that hook, which is faster when the material's own 3D tangent is expensive to form. It is a separate declaration because a generic fallback for compute_PK1_3D exists, so its presence alone cannot express the intent. A SmallStrain material ignores this and is always reduced from material_response.

compute_PK1_3D(material, F, dt, state)

Return the 3D first Piola-Kirchhoff stress for the dimensionality wrappers.

Extension point. A cross-package hook, optional for every material. The wrappers reach it only when the material declares tangent_source= compute_PK1_3D, and differentiate it for the reduced tangent. An AbstractHyperelastic subtype declares it and derives P from the strain energy Ψ.

update_state_from_3D!(state, material, F, dt)

Write a wrapped material's current/trial variables for the converged 3D deformation gradient F; for an immutable state, return its updated value.

Extension point. A cross-package hook, optional for every material model. The fallback takes the trial state from material_response using the declared strain measure, and performs no update for a stateless material. The history commit still happens through update_state! and update_states!.

deformation_gradient(cellvalues, qp, u_local)

Compute the deformation gradient F = I + ∇u at quadrature point qp from the given cellvalues and the element-local displacement vector u_local.

cellvalues is a Ferrite CellValues. qp is the quadrature point index. u_local is the displacement vector restricted to the current cell, in the order returned by celldofs(cell).

Returns a Tensor{2,dim} with the element type of u_local, where dim is the spatial dimension. For Float64 displacements this is Tensor{2,3,Float64,9} in 3D or Tensor{2,2,Float64,4} in 2D.

RemoteAssemblyFailure()

Recoverable assembly failure marker that try_stiffness_matrix returns on MPI ranks which did not fail locally, when another rank reported a LocalAssemblyFailure.

The rank where the failure occurred reports the original material-specific exception as result.error. All ranks receive converged = false, so an adaptive outer solver can reject the step collectively and retry with a smaller dt.

PlaneStressConvergenceError

Recoverable failure raised by PlaneStress when the local out-of-plane Newton solve or tangent condensation cannot produce a valid plane stress state. Catch this directly, or call try_stiffness_matrix to receive it as result.error.

VEPD_Detrez2010ConvergenceError

Recoverable failure raised when the local VEPD Detrez plastic correction does not converge within its fixed Newton iteration budget.

VEVP_MOAMMMConvergenceError

Recoverable failure raised when the local VEVP_MOAMMM plastic correction does not converge within its fixed Newton iteration budget.

alpha_value(alphavalues, qp::Int) -> Real

Return the scaling factor for quadrature point qp.

The element routine multiplies this factor into the integration weight before assembling into the residual and the stiffness, and compute_stresses applies the same factor to its output:

α_dΩ = alpha_value(alphavalues, qp) * getdetJdV(cellvalues, qp)

create_alpha_values creates alphavalues, and the assembler forwards it as the positional av argument to every material's element routine.

When create_assembler is called with ah=nothing, alphavalues is nothing and the built-in fallback returns 1.0, so no scaling is applied.

The alpha_value name, the AlphaValues companion type that create_alpha_values resolves by name, and the α_dΩ integration weight come from the CapriccioSimulation FE–MD coupling code: J. Roksvaag. CAPRICCIO - Tool to run concurrent Finite Element-Molecular Dynamics Simulations (Version 3.0.0), Zenodo (2026). https://doi.org/10.5281/zenodo.18326736

create_alpha_values(ah, cellvalues)

Build the per-cell alpha-values object that the assembler passes to alpha_value at every quadrature point.

The returned object av must satisfy alpha_value(av, qp::Int) -> Real. The assembler calls create_alpha_values once per cell during construction and stores the result for reuse in every assembly pass.

How to extend

Provide one of the following; the first match wins:

1. A method FerriteSolidMechanics.create_alpha_values(ah::MyAlphaSource, cellvalues)
2. A constructor AlphaValues(ah, cellvalues) in the same module that defines typeof(ah), which the default fallback calls automatically

create_alpha_values throws an ArgumentError if it finds neither.

When create_assembler is called with ah=nothing, create_alpha_values receives nothing and returns nothing. That pairs with the alpha_value(::Nothing, ::Int) = 1.0 fallback, so no scaling is applied.