Adaptive refinement: error estimation and marking

This page documents the theory behind marking for the dynamic hanging-node \(h\)-adaptivity modifier (hadaptivity): how EdelweissFE decides which elements to refine at a topology check, with particular attention to gradient-enhanced (implicit-gradient) damage models. The mechanics of how an element is then subdivided – octree children, curved-edge honouring, state transfer, and the exact hanging-node multi-point constraint – are documented under Model modifiers; this page is only about the refinement criterion.

Every formula stated here is implemented literally in marking; the recovered-gradient error indicator has been verified against a reference implementation and against analytic fields (exactly-linear field \(\to\) zero indicator; symmetric mesh and field \(\to\) equal indicators on mirror elements).

The adaptive loop and where markers sit

When the topology update runs is the solver’s decision, and the two solver families differ: the implicit solvers run update() at the start of every increment, while the explicit dynamic solver runs it at the end of a finalized increment and only every topology-check-frequency increments – see Adaptivity under explicit dynamics below. Either way, the update asks the modifier to plan() and then to apply() its decision. That update

1. evaluates each configured marker, obtaining from each a set of elements to refine,

2. unions those sets (see Reactive and predictive marking: why more than one marker),

3. keeps only elements below maxLevel,

4. subdivides them, rebalances, transfers state, and rebuilds the equation system.

A marker is a small object with a single responsibility – given the current model, return the set of elements it wants refined. The base class is MarkerBase; a marker is declared with a >>marker sub-keyword of *modelModifier, type=hAdaptivity. The available types are:

elementSet / nodeSet / surface

Purely geometric: mark a named set of elements, every element touching a named node set, or every element carrying a named surface. Typically used initialOnly for a fixed pre-refinement.

fieldOutput

Mark by thresholding an already-declared perElement *fieldOutput: operator (one of >, >=, <, <=, ==, !=) against threshold, applied element-wise, marking the element if any entry passes. The quantity being thresholded is shaped by that fieldOutput’s own f(x) – e.g. f(x)='np.abs(x)' for a magnitude, or f(x)='eigVal(x.reshape(-1,6))[:,0]...' for a principal stress – so the reduction lives in the (reusable) fieldOutput and the marker carries only the refinement policy. See FieldOutputMarker.

recoveryError

A recovery-based a posteriori error estimator on the gradient of a nodal field, with Dörfler bulk marking – the subject of most of this page. See RecoveryErrorMarker.

Adaptivity under explicit dynamics

The explicit dynamic solver NED (Solvers) also runs modifiers, but neither every increment nor at the same point in one, for reasons that are worth stating because they constrain how a marker may be configured.

Cadence. An explicit increment is thousands of times cheaper than an implicit one and there are correspondingly more of them – millions in a routine analysis. Evaluating a marker in each is not affordable, so the check runs every topology-check-frequency increments, and topology-check-frequency must be a multiple of output-frequency: a marker reads the last finalized field output, so a check on an increment whose output was never finalized would decide on stale results. Setting it to 0 runs the topology update exactly once, before the increment loop, which is what an initialOnly pre-refinement needs; a modifier that acts later than simulation start is then refused rather than silently ignored (see actsOnlyAtSimulationStart).

Placement. The check sits at the end of a finalized increment, not the start, for three independent reasons: the marker refines on the last finalized field output, so anywhere earlier it would decide on stale results; the velocity of a central-difference scheme is staggered half an increment behind the displacement, so the pairing of the two is unambiguous only between increments; and it keeps a topology change disjoint from the failure path, whose restore vectors are sized by the old equation system. Increment 0 is excluded deliberately – the time stepper yields a zero-length increment first, before anything has been solved, and a marker evaluated there would refine on the initial condition.

What a refinement changes for an explicit run. Refining lowers the stable time increment, so the critical time step is recomputed from the new mesh (unlike a mere contact-connectivity change, which leaves the mesh alone and therefore reuses it). The lumped mass is rebuilt; because that mass is the operator an explicit increment divides by, the solver reports the total mass, the per-component linear momentum and the kinetic energy across every topology change, and accumulates the relative mass drift over all of them so that many individually-tolerable changes cannot silently add up to a meaningful one. The velocity field is carried across refinement by the same warm-start interpolation as the displacement – omitting it would restart the new nodes from rest.

Motivation: gradient-enhanced damage

In an implicit-gradient (Peerlings/Geers) damage model the local driving field \(e_\mathrm{loc}\) (an equivalent strain, a plastic hardening variable, an energy release rate) is regularised into a nonlocal field \(\bar e\) by a screened-Poisson (Helmholtz) equation with internal length \(\ell\),

(1)\[\bar e - \ell^2\,\nabla^2 \bar e = e_\mathrm{loc} \qquad\text{in }\Omega,\qquad \nabla\bar e\cdot\mathbf n = 0 \text{ on }\partial\Omega,\]

and damage \(\omega\) evolves from the history \(\bar\kappa = \max_\tau \bar e(\tau)\). The nonlocal field \(\bar e\) is smooth (the Helmholtz operator is smoothing), but across the process zone it develops a steep gradient over a band of width \(\sim\ell\). Resolving that band – not the value of \(\bar e\), but its gradient – is what controls mesh objectivity of the dissipated energy.

Equation (1) is, term for term, a stationary heat-conduction problem with a reaction term: the “flux” conjugate to \(\bar e\) is \(\mathbf q = \ell^2\nabla\bar e\). Every error estimator developed for the temperature field of a thermal analysis therefore transfers directly by recovering \(\nabla\bar e\) in place of the heat flux – which is exactly what the recoveryError marker does. In EdelweissFE the nonlocal field is an ordinary nodal field named nonlocal damage, carried (for GC3D20/GC3D20R) by all 20 serendipity nodes at equal order with the displacement field.

Recovery-based a posteriori error estimation

For a finite-element field \(\bar e^h\) the discrete gradient \(\nabla\bar e^h\) is discontinuous across element boundaries and, for a smooth exact solution, converges one order more slowly than \(\bar e^h\) itself. The Zienkiewicz–Zhu idea is to construct a recovered gradient \(\nabla\bar e^{*}\) that is continuous and a better (super-convergent) approximation of the true gradient, and to use the discrepancy between the two as a local error indicator:

(2)\[\eta_K = \bigl\lVert\, \nabla\bar e^{*} - \nabla\bar e^{h}\, \bigr\rVert_{L^2(K)} = \left( \int_K \lvert \nabla\bar e^{*} - \nabla\bar e^{h}\rvert^2 \,\mathrm dV \right)^{1/2}, \qquad \eta^2 = \sum_K \eta_K^2 .\]

Because \(\nabla\bar e^{*}\) is closer to the exact gradient than \(\nabla\bar e^{h}\), the element quantity \(\eta_K\) estimates the local gradient (energy-norm) error. The intuition is purely geometric:

• where the mesh resolves the field, \(\nabla\bar e^{h}\) is already nearly continuous and close to \(\nabla\bar e^{*}\), so \(\eta_K\) is small;

• where a steep band sits on a too-coarse mesh, \(\nabla\bar e^{h}\) jumps sharply from one element to the next, the recovered (smoothed) field departs strongly from it, and \(\eta_K\) is large.

The constant factor \(\ell^2\) relating \(\nabla\bar e\) to the physical flux is omitted: it scales every \(\eta_K\) equally and so cannot change the ranking that marking is based on.

The integral in (2) is evaluated at the element’s \(2\times2\times2\) Gauss points, and both fields are sampled there: \(\nabla\bar e^{h}\) directly from the shape-function derivatives, and \(\nabla\bar e^{*}\) by interpolating the recovered nodal gradients with the element shape functions.

Gradient recovery

Two recovery schemes are available, selected by the recovery option.

Nodal averaging (recovery='averaging', Zienkiewicz–Zhu 1987). The recovered gradient at a node is the volume-weighted mean of the FE gradient sampled at that node by every element sharing it,

\[\nabla\bar e^{*}(\mathbf x_a) = \frac{\sum_{K \ni a} V_K\, \nabla\bar e^{h}\big|_K(\mathbf x_a)} {\sum_{K \ni a} V_K},\]

and \(\nabla\bar e^{*}\) elsewhere is the shape-function interpolation of these nodal values. Cheap and robust; the default. See _recover_averaging().

Superconvergent patch recovery (recovery='spr', Zienkiewicz–Zhu 1992). For each corner (vertex) node \(a\) a patch \(\omega_a\) is formed from the elements sharing it, and a complete polynomial of order \(p\) (matching the element order) is fitted to the FE gradient sampled at the super-convergent points of the patch – the element Gauss points – by least squares:

(3)\[\nabla\bar e^{h}(\mathbf x) \approx \mathbf p(\mathbf x)^\mathsf T\,\mathbf a, \qquad \mathbf a = \Bigl(\sum_{i\in\omega_a} \mathbf p_i\,\mathbf p_i^\mathsf T\Bigr)^{-1} \sum_{i\in\omega_a} \mathbf p_i\, \nabla\bar e^{h}_i, \qquad \mathbf p = [\,1,\,x,\,y,\,z,\,x^2,\,y^2,\,z^2,\,xy,\,yz,\,zx\,]^\mathsf T .\]

The recovered value at the node is the fit evaluated there, \(\nabla\bar e^{*}(\mathbf x_a) = \mathbf p(\mathbf x_a)^\mathsf T\mathbf a\) (with the fit centred on the node, simply the constant term \(\mathbf a_0\)). Edge-midside nodes, which are not patch centres, are recovered by averaging the fits of the two corner patches at each end of the edge, evaluated at the midside location. Boundary patches with too few sampling points for a full quadratic fall back to a linear fit, and, failing that, to the plain mean. See _recover_spr().

SPR is the sharper choice for the serendipity hexahedron and typically has near-unit effectivity (estimated / true error \(\approx 1\)), at the cost of a small least-squares solve per patch.

Note

Why the Gauss points. For the 20-node serendipity hexahedron the gradient is super-convergent at the \(2\times2\times2\) reduced-integration (Barlow) points – exactly the points at which GC3D20R integrates. Sampling the FE gradient there, rather than at the nodes, is what makes the recovered field higher-order. Both recovery schemes and the error integral use these points.

Marking: which of the ranked elements to refine

The indicators \(\{\eta_K\}\) are turned into a marked set by Dörfler (bulk) marking: sort the elements by \(\eta_K\) in descending order and mark the smallest set \(\mathcal M\) whose squared indicators first reach a fraction \(\theta\) of the total,

(4)\[\sum_{K\in\mathcal M}\eta_K^2 \;\ge\; \theta \sum_{K}\eta_K^2, \qquad \theta = \texttt{markFraction}\in(0,1].\]

This concentrates refinement on the largest error contributors. The marked set is additionally hard-capped at a fraction maxRefinedFraction of the eligible elements.

Why a plain ranking suffices here. Classical adaptivity chooses both where and how deeply to refine, iterating until an error target is met. The hanging-node modifier is normally run at a fixed single level (maxLevel=1) and a fixed subdivision factor (splitFactor): each marked element is split once, into \(\texttt{splitFactor}^3\) children. With the refinement depth fixed, only the location is adaptive – so marking collapses to “rank the elements and refine the worst ones”, which is precisely (4). A single-level split with splitFactor=3 (a \(3\times3\times3 = 27\)-child subdivision, exactly hanging-node coupled) delivers a three-fold local resolution increase in one step, often all that is affordable.

Why the budget cap matters with a direct solver. EdelweissFE assembles a single sparse system and factorises it with a direct solver (PARDISO). Each refinement changes the sparsity pattern, forcing a fresh symbolic factorisation, and direct-solve cost grows super-linearly in the number of unknowns. maxRefinedFraction bounds how many elements a single pass may add, keeping the per-increment factorisation cost predictable.

Reactive and predictive marking: why more than one marker

Markers are composable: one topology update refines the union of all configured markers’ sets. For a stationary problem one recoveryError marker is enough. For a propagating localization – the usual case in damage – a single marker is not, and the recommended configuration uses two markers together, for a reason worth stating carefully.

The recovery estimator is reactive: \(\eta_K\) only becomes large once an element already contains a steep, under-resolved gradient. But the process zone moves. By the time \(\eta_K\) flags an element, the damage band has already entered it on the coarse mesh – and a coarse discretization of a localizing band biases the dissipated energy, so the very mesh objectivity that gradient enhancement is meant to restore is compromised for that increment, and the band can “lock” onto the coarse discretization.

The remedy is a predictive marker that refines ahead of the front, keyed on a quantity that becomes active slightly before the nonlocal gradient sharpens – the local driving field \(e_\mathrm{loc}\), a plastic/damage hardening variable (e.g. alphaP), or a stress/strain criterion – expressed as a fieldOutput threshold marker. The two markers answer complementary questions:

Marker	Question it answers	Type
Reactive	Where is the solution under-resolved now?	recoveryError (this page)
Predictive	Where is the solution about to develop a feature?	fieldOutput threshold on the driving field

Neither alone suffices. Reactive only is always one step behind the front: the band localizes on a coarse element before refinement catches up, and, with a direct solver, catching up costs extra corrector re-solves (each a fresh factorisation). Predictive only is a crude heuristic – no error control, it over-refines wherever damage merely initiates, cannot adapt to the actual discretization error, and is blind to any feature not tied to that one state variable. Their union gives both: the predictive term keeps the mesh fine before the band arrives (preserving mesh objectivity and avoiding extra factorisations), and the reactive estimator supplies principled, error-controlled refinement of the active zone and anywhere else the discretization is poor. This is the predictor–corrector strategy of phase-field fracture adaptivity, expressed as two composable markers.

Usage

A recoveryError marker names the nodal field whose recovered-gradient error drives refinement, and optionally the recovery method, bulk fraction and budget cap:

*modelModifier, type=hAdaptivity, name=amr
** reactive: ZZ recovered-gradient error of the nonlocal field localizes on the process zone
>>marker, type=recoveryError, nodeField='nonlocal damage', markFraction=0.5, maxRefinedFraction=0.1, recovery=spr
** predictive: refine ahead of the front where the local driving field starts to grow
>>marker, type=fieldOutput, fieldOutput=alphaPForAMR, operator='>', threshold=1e-4
elSet=concrete
maxLevel=1
splitFactor=3

The predictive marker’s fieldOutput (here alphaPForAMR) is an ordinary perElement output of the local driving variable, covering the quadrature points of interest with no f(x) reduction. Both markers accept a halo option (0 by default): a number of node-adjacent-element rings added around the marked set, bounded to the refineable elements – useful when a bare threshold would leave a conforming band’s immediate surroundings coarse, e.g. >>marker, type=fieldOutput, fieldOutput=alphaPForAMR, operator='>', threshold=1e-4, halo=1. The full option set of the recoveryError sub-keyword:

``hAdaptivity`` options (the ``>>marker`` sub-keyword is rendered as its own table):

Option	Type / Default	Description
moduleOptions	dict, default={}	Internal
elSet	str, default=None	Fallback for ‘refineElSet’ if that is not given. Does not restrict marking itself – each ‘>>marker’ scopes its own eligible elements (a fieldOutput’s associated set, an elementSet/nodeSet/surface’s members).
refineElSet	str, default=None	Restrict the AMR octree mirror itself to this element set, e.g. the solid elements in a mesh that also contains contact-facet elements. Elements outside this set never become octree roots and are left untouched by refinement. Defaults to ‘elSet’ if given, otherwise to every 20-node (HEX20-family) element in the model.
maxLevel	int, default=1	Maximum refinement level.
minMarkedElements	int, default=1	Minimum number of eligible elements one evaluation of the markers must mark for a refinement pass to happen; with fewer, nothing is refined and no equation system is rebuilt. Marks are not carried over to later evaluations: they are a function of the current state only. Individual markers may cap their own marks (e.g. ‘maxRefinedFraction’) or expand them (e.g. ‘halo’) before this count is taken. Default 1 refines as soon as any element is marked.
splitFactor	int, default=2	Number of equal parts per axis a marked element is split into (2 = octree bisection into 8 children; 3 = 3x3x3 = 27 children, etc.). The hanging-node coupling stays exact for any factor.
elementType	str, default=None	Element type to instantiate for children (default: like parents).
elementProvider	str, default=’marmot’	Element provider.
stateTransfer	str, default=’nearestQp’	Quadrature-point state-transfer strategy for the whole state block: nearestQp|projection|virgin.
stateTransferOverrides	str, default=None	Per-state-variable overrides routing named variables to a different strategy, e.g. ‘strain:projection, stress:virgin’. Comma-separated ‘name:strategy’ pairs.
marker	tuple, default=()	AMR marker definition. At least one is required.

Keyword: marker

Option	Type / Default	Description
type	str (required)	Marker type, resolved through the ‘marker’ registry: fieldOutput, elementSet, nodeSet, surface, recoveryError. The type-specific options are defined by the selected marker’s own schema.
initialOnly	bool, default=False	Evaluate only once at simulation start

A minimal, self-contained regression case (a two-element GC3D20R cantilever whose nonlocal-damage gradient concentrates at the fixed end, so the marker refines exactly one element) is:

testfiles/marmot/AMR_RecoveryError/test.inp (recovery='averaging'; the AMR_RecoveryErrorSPR sibling is identical but for recovery=spr).

** Dynamic h-adaptivity regression test for the recoveryError marker: a Zienkiewicz-Zhu recovered-
** gradient error indicator on the nonlocal damage field, with Doerfler bulk marking. Two coarse
** GC3D20 elements form a cantilever; the outer element is stretched 5% (x in [1, 2.05] instead of
** [1, 2]) to break an exact tie in the two elements' standard-averaging ZZ error under this
** material/loading -- with an unstretched (exactly symmetric) mesh both elements recover an
** identical indicator, so which one wins the marking decision is an environment-dependent
** argsort tie-break rather than physics. With the stretch, bending concentrates the gradient of
** the nonlocal damage field near the fixed end as intended, so the recovery marker (worst element
** by ZZ error, capped by the DOF budget) unambiguously refines exactly one element -- creating a
** 2:1 hanging interface the MPC couples across both fields (displacement + nonlocal damage). Pins
** the dynamic solution against U.ref, and thus exercises the full modelModifier path: recovery
** marking -> octree refinement -> model mutation -> warm start -> quadrature-point state transfer
** -> hanging-node MPC -> rebuild -> solve.
*material, name=AT2PHASEFIELD, id=mat
20000, 0.2, 0.1, 5

*section, name=sec, material=mat, type=solid
all

*modelModifier, type=hAdaptivity, name=amr
>>marker, type=recoveryError, nodeField='nonlocal damage', markFraction=0.5, maxRefinedFraction=0.5
maxLevel=1

*job, name=amr, domain=3d
*solver, solver=NIST, name=S
*fieldOutput
>>perNode, elSet=all, field=displacement, result=U, name=U
>>perNode, elSet=all, field=nonlocal damage, result=U, name=D

*step, solver=S
maxInc=0.5, minInc=1e-4, maxNumInc=50, maxIter=25,
>>dirichlet, name=x0, nSet=x0, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=x1, nSet=x1, field=displacement, 3=0.0005

*Node
1, 0.0, 1.0, 1.0
2, 0.0, 0.0, 1.0
3, 0.0, 0.0, 0.0
4, 0.0, 1.0, 0.0
5, 1.0, 1.0, 1.0
6, 1.0, 0.0, 1.0
7, 1.0, 0.0, 0.0
8, 1.0, 1.0, 0.0
9, 0.0, 0.5, 1.0
10, 0.0, 0.0, 0.5
11, 0.0, 0.5, 0.0
12, 0.0, 1.0, 0.5
13, 1.0, 0.5, 1.0
14, 1.0, 0.0, 0.5
15, 1.0, 0.5, 0.0
16, 1.0, 1.0, 0.5
17, 0.5, 1.0, 1.0
18, 0.5, 0.0, 1.0
19, 0.5, 0.0, 0.0
20, 0.5, 1.0, 0.0
21, 2.0500, 1.0, 1.0
22, 2.0500, 0.0, 1.0
23, 2.0500, 0.0, 0.0
24, 2.0500, 1.0, 0.0
25, 2.0500, 0.5, 1.0
26, 2.0500, 0.0, 0.5
27, 2.0500, 0.5, 0.0
28, 2.0500, 1.0, 0.5
29, 1.5250, 1.0, 1.0
30, 1.5250, 0.0, 1.0
31, 1.5250, 0.0, 0.0
32, 1.5250, 1.0, 0.0
*Element, type=GC3D20
1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
2, 5, 6, 7, 8, 21, 22, 23, 24, 13, 14, 15, 16, 25, 26, 27, 28, 29, 30, 31, 32
*elSet, elSet=all, generate
1, 2, 1
*nSet, nSet=x0
1, 2, 3, 4, 9, 10, 11, 12
*nSet, nSet=x1
21, 22, 23, 24, 25, 26, 27, 28
*nSet, nSet=yf
2, 3, 6, 7, 10, 14, 18, 19
22, 23, 26, 30, 31
*nSet, nSet=zf
3, 4, 7, 8, 11, 15, 19, 20
23, 24, 27, 31, 32

Implementation

Refinement marking (WS-G): the markers that decide which elements an adaptivity mechanism refines.

Each marker is a small, registry-resolved policy object (see edelweissfe.config.markerlibrary) implementing MarkerBase.mark(). Topological markers select from a set/surface; the FieldOutputMarker thresholds a per-element quantity. Deliberately, a marker never re-derives a field quantity itself: the reduction (a magnitude, a principal stress via eigVal, a norm, …) lives in the referenced fieldOutput’s own f(x) – one first-class, reusable quantity – and the marker only applies the refinement decision (a threshold) on top of it.

class edelweissfe.adaptivity.marking.FieldOutputMarker(fieldOutputName, threshold, operator='>=', halo=0, initialOnly=False)

Marks every element whose fieldOutput result satisfies value <operator> threshold at any entry (any quadrature point / component).

The marker is a pure threshold decision: the quantity it compares is produced by the referenced perElement fieldOutput – either raw, or shaped by that fieldOutput’s own f(x) (abs(x) for a magnitude, eigVal for a principal stress, a norm, …). Keeping the reduction in the fieldOutput (where it is a first-class, reusable/visualizable quantity) and only the threshold policy here is the division of labour – the fieldOutput says what the quantity is, the marker says when it is large enough to refine. Marking is driven by the exact numbers the fieldOutput reports, and an element is marked if np.any of the comparison over its result row holds, so a multi-QP / multi-component row is decided by its worst entry.

Parameters:

halo – Number of rings of node-adjacent elements to add around the thresholded set (see _growByNeighbors()), the same buffer mechanism as RecoveryErrorMarker’s own halo. Applied after the threshold decision, so the final marked count can exceed what the bare threshold would have produced. Grows only within the refineable elements, so it never leaks onto elements outside refineElSet. 0 (default) keeps the bare thresholded set.

schema

alias of FieldOutputMarkerSchema

classmethod fromOptions(options)

Construct this marker from its raw parsed >>marker options.

options is the >>marker block’s option mapping with the type dispatch key already removed (the adaptivity mechanism uses type to look the class up before calling this). The options are validated and coerced against schema via buildSchemaFromOptions(), then mapped onto the constructor.

checkReferences(model)

Raise if a set, surface or fieldOutput this marker reads does not exist in model.

For a caller that knows all of them must exist when the run starts (e.g. the Schur complement of a composed solver, see edelweissfe.linsolve.composition.schurcomplement), so that a typo fails at setup instead of silently marking nothing. Nothing to check by default.

Parameters:

model – The FEModel.

isReady(model)

Whether the marker has what it decides on, e.g. a fieldOutput result, which exists only after the first converged increment. Always, by default.

Parameters:

model – The FEModel.

Return type:

bool

mark(model, candidateElements, mesh)

Select elements: those to refine, for adaptive mesh refinement.

Parameters:

• model – The FEModel, for the sets, surfaces, fields and fieldOutputs a marker reads.

• candidateElements – The elements the marker may select (e.g. the refineable ones), with their node adjacency.

• mesh – The adaptivity mechanism’s mesh mirror; None outside adaptive refinement (no marker here reads it).

Returns:

The selected elements.

Return type:

set

selectionStamp(model)

A value that changes whenever what isReady() and mark() read can change, other than the model itself (its elements, sets and surfaces, which change only with a new model): STABLE_FOR_THE_MODEL for a marker that reads nothing else, or None if the marker cannot tell. A caller may keep a selection while the stamp stays the same (the Schur complement of a composed solver does); None, the default, means: mark again.

Parameters:

model – The FEModel.

class edelweissfe.adaptivity.marking.RecoveryErrorMarker(nodeFieldName, markFraction=0.5, maxRefinedFraction=0.1, fieldThreshold=1e-06, relTol=0.001, halo=0, recovery='averaging', entry='U', initialOnly=False)

Zienkiewicz-Zhu recovery-based error indicator on the gradient of a nodal field, with Doerfler bulk marking. Aimed at gradient-enhanced damage: the nonlocal driving field \(\bar\varepsilon\) (governed by a screened-Poisson/Helmholtz equation) is smooth, but its gradient localizes in the process zone, and that is where the mesh must be resolved.

The per-element indicator is the L2 norm of the difference between a recovered, continuous gradient field and the raw (discontinuous) FE gradient,

\[\eta_K = \lVert \nabla\bar\varepsilon^{*} - \nabla\bar\varepsilon^{h} \rVert_{L^2(K)},\]

which is exactly the thermal-flux ZZ estimator applied to the nonlocal field (the constant length-scale factor \(\ell^2\) cancels out of a pure ranking, so it is omitted). The recovered field \(\nabla\bar\varepsilon^{*}\) is built by volume-weighted nodal averaging of the FE gradient sampled at the 2x2x2 (superconvergent) Gauss points – the serendipity HEX20 sweet spot. (An SPR patch-least-squares variant is a planned upgrade behind recovery.)

Selection is by Doerfler bulk marking on the squared indicators: the worst elements whose \(\sum \eta_K^2\) first reaches markFraction of the total are marked – but never more than maxRefinedFraction of all elements. Under a fixed single-level refinement this is the natural criterion: refinement depth is fixed, so marking reduces to ranking, and the hard cap keeps the direct-solver (PARDISO) factorization cost bounded.

Doerfler is a purely relative criterion, so on its own it always marks a markFraction share of the ranking – even in the elastic phase, before any damage, where the nonlocal field is machine zero (~1e-30). Crucially, a relative gate cannot fix this: the noise-level field’s recovered-vs-FE gradient mismatch is O(1) relative to itself, indistinguishable from a genuinely localizing field. Marking is therefore gated by an absolute floor fieldThreshold on the peak nodal |field| per element – an element is a candidate only where the driving field is significant (the active-damage field is 1e-5 and up, ~25 orders above the elastic noise, so the separation is unambiguous). Among those significant elements a secondary ZZ global relative error gate relTol skips regions the mesh already resolves well.

Parameters:

• nodeFieldName – Name of the nodal field whose recovered-gradient error drives marking (e.g. 'nonlocal damage').

• markFraction – Doerfler bulk fraction theta in (0, 1]: mark the worst elements accumulating this fraction of the total squared error.

• maxRefinedFraction – Hard cap on the fraction of eligible elements a single pass may mark.

• fieldThreshold – Absolute floor on the peak nodal |field| per element: an element is a refinement candidate only where the driving field itself is significant. This is what keeps the elastic phase (where the nonlocal field is machine zero, ~1e-30) from being refined – a purely relative criterion cannot, because the noise-level field’s recovered-vs-FE gradient mismatch is O(1) relative to itself. Set it to a small fraction of the field value at damage initiation (for a strain-like driving field, damage onset is O(f_t / E), so a default of 1e-6 sits well below onset yet ~20 orders above the elastic noise). Set to 0.0 to disable and refine on gradient roughness alone.

• relTol – Global relative-error threshold, applied among the significant elements only: while \(\lVert \eta \rVert / \lVert \nabla\bar\varepsilon^{h} \rVert < \text{relTol}\) the field is treated as well resolved and no element is refined. Set to 0.0 to disable.

• halo – Number of rings of node-adjacent elements to add around the Doerfler-marked set (a refinement buffer). The ZZ estimator responds to gradient error, which peaks on the flanks of a localization band and dips over its (locally smooth) core, so on its own it can refine two rows straddling a crack while leaving the core coarse. halo=1 bridges those rows across the core; a larger halo also refines ahead of the propagating front. 0 (default) keeps the bare Doerfler set. The halo grows only within the refineable elements, so it never leaks elsewhere, and it is applied after the maxRefinedFraction cap, so the final count may exceed it.

• recovery – Recovery method: 'averaging' (volume-weighted nodal averaging, ZZ 1987) or 'spr' (superconvergent patch recovery, ZZ 1992 – the sharper choice for serendipity elements).

• entry – Node-field value entry to read, one of ENTRIES ('U', the current converged field, by default).

schema

alias of RecoveryErrorMarkerSchema

classmethod fromOptions(options)

Construct this marker from its raw parsed >>marker options.

options is the >>marker block’s option mapping with the type dispatch key already removed (the adaptivity mechanism uses type to look the class up before calling this). The options are validated and coerced against schema via buildSchemaFromOptions(), then mapped onto the constructor.

mark(model, candidateElements, mesh)

Select elements: those to refine, for adaptive mesh refinement.

Parameters:

• model – The FEModel, for the sets, surfaces, fields and fieldOutputs a marker reads.

• candidateElements – The elements the marker may select (e.g. the refineable ones), with their node adjacency.

• mesh – The adaptivity mechanism’s mesh mirror; None outside adaptive refinement (no marker here reads it).

Returns:

The selected elements.

Return type:

set

Per increment the marker gathers, for every eligible 20-node element, its node coordinates and nonlocal-field values plus a compact global node index, then evaluates the indicator (2) for the whole element set at once in _recovery_indicators(): all per-element \(3\times3\) Jacobians, gradients and the nodal-averaging accumulation are batched into single numpy calls (einsum and a stacked linalg.inv), so there is no Python-level per-element loop for the default path. The recovery step is delegated to _recover_averaging() or _recover_spr(), and the Dörfler selection (4) and budget cap are applied to the resulting ranking.

Limitations and extensions

• The recovered field is resolved at the fixed single refinement level; the estimator ranks where to refine but does not itself choose refinement depth (that is maxLevel / splitFactor).

• SPR uses a complete quadratic patch basis, reduced to linear (or a mean) on patches too small for it; an anisotropic (Hessian-metric) recovered field, which would allow band-aligned stretched refinement, is not implemented – the modifier refines isotropically.

• Refinement is monotone: there is no de-refinement (coarsening) of the wake behind a passed front.

• The internal length \(\ell\) sets the resolution the band actually needs (roughly \(h \lesssim \ell/2\) for the quadratic serendipity element); the marker chooses where to reach that resolution, but the initial mesh must at least marginally resolve \(\ell\) for the estimator to fire on a band that is otherwise never captured.

References

• O. C. Zienkiewicz, J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Numer. Methods Eng. 24 (1987) 337–357.

• O. C. Zienkiewicz, J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique, Int. J. Numer. Methods Eng. 33 (1992) 1331–1364.

• W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996) 1106–1124.

• R. H. J. Peerlings, R. de Borst, W. A. M. Brekelmans, M. G. D. Geers, Gradient-enhanced damage modelling of concrete fracture, Mech. Cohes.-Frict. Mater. 3 (1998) 323–342.

• T. Heister, M. F. Wheeler, T. Wick, A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach, Comput. Methods Appl. Mech. Eng. 290 (2015) 466–495.