The topology pipeline: how a model changes itself

An adaptive analysis mutates its own model while it runs. Refinement splits elements, a deposition module adds them, contact and tie surfaces regenerate their facets. Every one of those creates and destroys elements and nodes – and every one of them must produce identical results when the analysis is resumed from a checkpoint, or restored material state lands on the wrong element.

This page describes how EdelweissFE arranges that. It is the page to read before writing a new model modifier, or before writing anything that caches data derived from the mesh.

The short version, and the mnemonic worth remembering:

Modifiers plan and apply; mesh-dependents refresh; then we solve.

One increment, three phases

╔══════════════════════════════════════════════════════════[ window OPEN ]══╗
║ PHASE 1 — TOPOLOGY UPDATE   the only place elements/nodes are born or die ║
║                                                                           ║
║   amr ─────────▶ printer ─────────▶ facets ─────▶ any plan applied?       ║
║   plan(model)    plan(model)        plan(model)          │   │            ║
║      ↓              ↓                  ↓                 │   │            ║
║   apply(…)       apply(…)           apply(…)             │   │            ║
║        modifiers run in declared order                   │   │            ║
║   ┌──────────────── yes: run another round ◀─────────────┘   │            ║
║   └──▶ (back to amr)                      no: fixed point ───┼──┐          ║
╚══════════════════════════════════════════════════════════════╪══╪═════════╝
                                                               │  │
╔══════════════════════════════════════════════════════════════▼══▼═════════╗
║ PHASE 2 — REFRESH  [ window CLOSED ]  pure readers: nothing born or dies  ║
║   each mesh-dependent asks: has topology.version moved since I last looked?║
║        yes → refresh once, from the NET change across all rounds           ║
║   order is irrelevant — none of them can invalidate another's work         ║
╚═══════════════════════════════════════════════════════════════════════════╝
                                      │
╔══════════════════════════════════════▼════════════════════════════════════╗
║ PHASE 3 — SOLVE   assemble the equation system and iterate to convergence  ║
╚═══════════════════════════════════════════════════════════════════════════╝

Phase 1 is update(), phase 2 is refreshMeshDependents(). The solver calls both, in that order, before every solve.

Why rounds

Modifiers depend on each other. Refinement invalidates a tied surface’s facets; a deposition modifier creates elements refinement may then want to split; a 2:1 balance may need to refine what another modifier just activated. A declared dependency order cannot express mutual dependence, and asking users to get one right is a poor trade.

So each round offers every modifier the net change since that modifier last planned. A round in which nobody plans anything is the fixed point, and it ends the phase. Determinism comes from the structure, not from luck: within a round, modifiers run in the order they were declared in the input file.

A modifier that keeps planning in response to its own output would never settle. That is a bug in the modifier, and maxRounds turns it into a loud error naming the offender rather than a hang.

Two modifiers must not change the same element within one update, in any round – whichever ran second silently wins, and which that is depends on declaration order. The pipeline refuses it, naming both modifiers, the element, and the rounds they touched it in. Reacting to another modifier’s change in a later round is fine; mutating the same element it mutated is not.

An apply that returns an empty ModelChange counts as nothing having happened: it is not recorded in the history, it does not bump the version, and it does not trigger an equation-system rebuild. Planning and then finding nothing left to do is a legitimate outcome, not a change.

Writing a model modifier: `plan` and `apply`

A modifier implements two halves, and the split is the single most important thing on this page.

plan(): Decide what to do, without doing it. May read anything: markers, node fields, the current time, the step. Returns a serializable description of the decision, or None if there is nothing to do.
apply(): Carry the decision out, mutating the model. Must not read solution state – it is a pure function of (model, plan). Everything the decision depended on belongs in the plan.

def plan(self, model, change, step):
    if change is not None and not change.touchesElementSet(self.myElSet):
        return None                      # not my business -- let the pipeline settle
    marked = self.evaluateMarkers(model) # reads solution state: allowed here
    return None if not marked else MyPlan(ids=sorted(marked))

def apply(self, model, plan):
    numbers = model.topology.reserveElementNumbers(len(plan.ids))
    for identifier, elNumber in zip(plan.ids, numbers):
        model.createElement(self.makeElement(elNumber, identifier))
    return change                        # a ModelChange describing what happened

Two rules follow from the signatures and are worth stating explicitly:

Return ``None`` from plan when the change does not touch your domain. That is what lets the pipeline reach a fixed point instead of looping. change is None on the first round of an update, meaning “evaluate freshly”.

Never write model.elements or model.nodes directly. Use reserveElementNumbers(), createElement() and removeElement(), and their node-side counterparts reserveNodeNumbers() and createNode().

Say so if your modifier is purely reactive. A modifier that can only act in response to another one’s mutation – the implicit surface-facet retiling is the only current example – sets initiatesTopologyChanges = False. A solver that never runs the topology update (the explicit ones) then still accepts a model carrying it, while still refusing one whose modifiers would silently never run.

Why the split exists: restart

A resumed run must rebuild exactly the mesh it resumed – same elements, same numbers, same connectivity – because element state is restored by number.

The obvious way to arrange that is for each modifier to serialize its own history and implement its own replay. EdelweissFE did exactly that once, and it produced a real bug: a resumed analysis rebuilt a differently-numbered mesh, so heavily damaged elements were restored virgin, and the solve collapsed into a cutback spiral several increments later.

The root cause was not the numbering. It was that replay ran through different code than the live run, and two implementations of one mutation always drift apart. The drift showed up as a handful of individually-reasonable patches – “skip the warm start on replay”, “this is never really the first call”, “materialise level-wise instead of occasion-wise” – which together changed the result.

The plan/apply split removes the second implementation:

    LIVE RUN                                      RESTART REPLAY
modifier.plan(model, …)                        read recorded plan
reads markers, state, time                     no markers re-evaluated
     │        │                                     │            │
     │        └── records ──▶ topology.history ─ restores ───────┘
     │                                                           │
     └──────────── applies ──▶┌───────────────────────┐◀─ applies┘
                              │ modifier.apply(model, │
                              │               plan)   │
                              │ ONE code path         │
                              │ reads no solution state│
                              └───────────────────────┘
                                          │
                    identical element numbers · identical topology

history records every applied decision – the modifier’s own serialized plan, via encodePlan(), plus the resulting fingerprint(). A restart replays it through replayHistory(), which calls the same apply. plan is never called during a replay.

The fingerprint is not decoration. Recorded per decision, it is what lets a resumed run prove it rebuilt the mesh it was checkpointed with. It is not free either: it walks the whole mesh, measured at 0.188 s on 64k elements / 69k nodes, and a live run pays that once per applied decision – not per iteration, but often enough to notice on a large model.

A replay pays it once: the recorded digests are carried forward into the replayed history, and the mesh is checked against the last one after the whole history has been applied. Checking every record instead would make replaying a long history O(records x mesh) rather than O(mesh) – a few hundred refinements on a mesh of tens of thousands of elements took minutes that way. When the final check does report a divergence, set verifyFingerprintsPerRecord and replay again: that turns “the resumed run diverged somewhere” into “it diverged at record 12, modifier amr , round 2” – a divergence you can bisect rather than hunt.

The fingerprint hashes every element’s number, type and connectivity and every node’s label and reference coordinates – the mesh, not the motion. The surface nodes of a discrete rigid body are the one exception to the rule that a node’s coordinates are its reference coordinates: the body writes its current position into them. It therefore reports their reference coordinates to the fingerprint (referenceCoordinatesOfMovedNodes()). Otherwise a decision recorded after the body had moved could never be verified on resume, because the replay runs before the displacement is restored, with the body still at its reference position.

For the same reason the node-field bookkeeping a mutator requests after each change – activating field variables on every node, resizing every NodeField, relinking every field variable – is deferred to the end of the replay window and run once (changes() with deferFieldBookkeeping). It is recomputed from the final mesh either way, so the result is the same; only the intermediate layouts, which no increment ever solves on, are skipped. The live per-increment window does not defer, and apply itself is unaware of the difference: it issues the same calls in both cases.

What a modifier carries between decisions is declared like every other component’s state, in checkpointedState – hAdaptivity’s flag that its step-start markers were evaluated, for instance – and restored after the replay. Everything else plan reads is the restored model and solution state. This is deliberately separate from apply: it affects how the next decision is made, so getting it wrong cannot corrupt the mesh.

Element and node numbers

Element numbers come from one monotonic allocator, reserveElementNumbers(). They are never recycled, and never derived from max(model.elements). Node labels get exactly the same treatment from reserveNodeNumbers(), for exactly the same reasons; the rest of this section says “element” for both.

That is not tidiness. Deriving the next number from the current maximum makes numbering a function of the deletion history as well as the creation history – a contact facet set that is deleted and rebuilt hands the freed numbers straight back out – so a replay would have to reproduce creations, deletions and their interleaving. With one monotonic counter it only has to reproduce the ordered creation sequence, which is exactly what the recorded history holds.

Never recycling also means a number refers to one element for the model’s entire lifetime, so a reference cached by number cannot silently come to mean something else.

Numbers may only be reserved inside a topology window (changes()), which the pipeline opens around phase 1. Outside it, creating or deleting an element or a node raises TopologyError. This is what makes “only model modifiers change the topology” an enforced property rather than a convention.

There is exactly one allocator per model, and everything that mints draws from it: the mesh generators (which run in a setup-time window opened around model setup), and the node registry of adaptive refinement. That registry lives in the model-agnostic octree layer, so it is handed reserveNodeNumbers() as a plain count -> range callable rather than a model. Without one it falls back to minting max + 1 itself, which keeps the octree usable standalone but is only safe while nothing else mints.

Reacting to a change: three tiers, deliberately

Anything that caches data derived from the mesh has to notice when the mesh changes. There are three ways, in increasing order of power – use the weakest one that works.

Tier 0 – do nothing. Node sets, element sets, surfaces and node fields have stable identity: they are mutated in place, never replaced. A component that caches model.nodeSets["mySet"] and iterates it later sees any refinement automatically. Most code needs nothing else, and that is by design.

Tier 1 – a set-changed self-check. A component that pre-sizes a derived array to a container’s size needs to notice a resize, since the array does not grow on its own. It compares the container’s version counter at its own point of use (ConstraintBase._checkSetChanged, StepActionBase._checkSetChanged). This has no timing dependency at all – it checks when it uses the data – which makes it correct on a replay for free.

Tier 2 – MeshDependent. For cached geometry: contact facets, tie records, projections. Register once with registerMeshDependent(), implement refresh(), and phase 2 calls it once per increment with the net change.

class MyConstraint(ConstraintBase, MeshDependent):
    def __init__(self, name, model, ...):
        model.topology.registerMeshDependent(self)

    def refresh(self, model, change):
        if not change.touchesSurface(self.mySurface):
            return False                 # not my business
        self.rebuildCachedGeometry(model)
        return True                      # my DOF footprint changed

Registration is the freshness guarantee, not ceremony: a consumer that is not registered is never told. That matters most for consumers the solver does not otherwise tick – multi-point constraints live in model.multiPointConstraints, which no per-increment sweep iterates, so a tie can learn about refinement only this way.

refresh must not create or delete entities. Topology is a modifier’s business, and the window is closed by the time phase 2 runs, so an attempt raises.

Why pull, and not a push notification

EdelweissFE once notified consumers synchronously, at the instant of each mutation. That mechanism is gone, and the reason is worth understanding, because it is a general one.

With modifiers running to a fixed point, a per-mutation callback necessarily fires mid-pipeline. Round 3 rebuilds what round 2 invalidated, so a consumer notified in round 1 is handed a state that no longer exists by the time the solve begins. And a consumer that mutates in response – every tie and contact does, since they mint facet elements – would do so re-entrantly, inside the modifier’s own loop.

Push is not merely riskier here; it is at the wrong granularity. It reports transients, and consumers want the settled model. So phase 2 pulls: each consumer diffs against its own last-seen version and reconciles once, from the change coalesced across every round.

A checklist for a new modifier

Implement plan (may read state) and apply (may not).
Return None from plan when the incoming change does not touch your domain.
Reserve element numbers from the model; never write model.elements.
Implement encodePlan/decodePlan so your decision survives a checkpoint.
Declare checkpointedState: what plan carries from one decision to the next (often nothing).
Verify with a restart round-trip and compare fingerprint(), not just element counts.