Second-order elements: a coarse mesh that beats a fine one

In a first-order (P1) triangle the vector potential is linear, so B — its gradient — is constant over the element. Every quantity an engineer actually cares about in a motor model is computed from B: air-gap flux density, torque, iron loss, back-EMF. A piecewise-constant B across a four-element-thick air gap is the quiet reason a torque number wobbles when you remesh.

The trade, stated honestly

A P2 triangle carries six nodes instead of three — the three corners plus one on the middle of each edge. A is quadratic, B is linear within the element, and the local stiffness block is 6×6 rather than 3×3. You pay roughly four times the unknowns per element and a denser matrix row. The question is whether you get more than four times the accuracy back, because if you do not, you should simply have refined the mesh instead.

The measurement that answers it

A 3-pole-pair, 18-slot PMSM in a nonlinear transient run with an anti-periodic sector, an arc sliding band and three coils. Four runs: coarse P1, coarse P2 (the same 6330-element mesh), and 4×-refined versions of both. Sampled on a 241-point circle in the air gap and in the iron, against the refined P2 reference:

Run	rel L2 in A	rel L2 in |B|, air gap
coarse P1	1.4e-3 … 2.6e-3	5.5e-3
4×-refined P1	0.9e-3 … 1.5e-3	5.2e-3
coarse P2	5e-5 … 2.5e-4	7e-5 … 6e-4

Coarse P2 is an order of magnitude closer to the converged answer than P1 on a mesh four times finer — and all four runs converged in the same number of Newton iterations. Notice also how little the refined P1 row moved: in the air gap, quadrupling the P1 mesh bought 5 %. That is the signature of an error that is not going to be refined away at any mesh size you would tolerate.

The same shape shows up in a clean analytic case — a round conductor in air with a circular Dirichlet boundary. L2(|B|) under refinement: P1 gives 1.38e-2 → 7.7e-3 → 3.6e-3; P2 gives 1.66e-3 → 6.6e-4 → 4.9e-4. P2 on the coarsest mesh beats P1 on the finest.

And on a quadratic patch test — a rectangle with Dirichlet and Neumann boundaries and a source whose exact solution is A = x² — P2 reproduces the exact answer to 1.9e-14 relative in A and 1.2e-15 in B, on every mesh, while P1 shows the expected O(h²) error. That is not an accuracy result so much as a correctness one: it is the test that catches a wrong shape-function gradient or a mis-numbered midside node.

Three details that decide whether you get the accuracy

Midside numbering, and gradients that are actually gradients

Which edge each midside node belongs to is a convention, and the two plausible conventions differ by a rotation. Get it wrong and the model still solves — it simply solves something else. A related failure mode is worse: an implementation that fakes the quadratic gradients by averaging corner values produces a 6×6 element matrix of rank 3, which is to say a P1 element wearing six nodes. Both are invisible without a patch test, which is why the A = x² case above matters more than it looks.

Material state is per quadrature point, not per element

A P1 element has one B, so it has one permeability. A P2 element has a B that varies across it, so the saturating material model has to be evaluated at each quadrature point. The subtle failure here is in the non-saturating path: if K is assembled once up front, with no Newton loop to trigger a field recovery first, every quadrature point still carries its default μ = 1 — and the stator and rotor cores behave as air. In our own case that collapsed peak |B| in a PMSM from 2.5 T to 0.99 T (the magnets alone). Saturating runs never showed it, because their Newton loop recovers the field before each re-assembly. A bug that only appears when you turn the nonlinearity off is exactly the kind a validation suite has to be built to catch.

Interfaces and post-processing have to be quadratic too

A P2 boundary ring alternates corner and midside nodes. Tying a periodic-boundary slave node linearly between two adjacent ring entries caps an otherwise third-order sector model at second order — all the cost of P2, some of the benefit. The same applies to a sliding band, which ties node values across a rotating interface: the moving ring has to be interpolated against the full quadratic master edge.

Post-processing is the one people forget. Integrating over the three corners only — the natural thing to do, since that is what P1 needs — injects an O(h²·∂²A) quadrature error, the same order as the P1 interpolation error, on top of a field accurate to about 1e-4. The quadrature then becomes the dominant error term and silently erases the accuracy you paid for. For a straight-sided quadratic triangle the corner shape functions integrate to zero and each midside function to S/3, so the correct element mean is (A₃+A₄+A₅)/3 over the midside values — the corners do not appear at all.

What P2 supports in Nabla

Planar and axisymmetric; static, transient and time-harmonic; linear and saturating; with eddy currents, with backward-Euler or Crank–Nicolson time integration, with circuit coupling (staggered and monolithic), with Dirichlet, Neumann and far-field boundaries, with periodic and anti-periodic interfaces, with rotary and linear motion through the averaging sliding band, and through the thermal module.

Two configurations are rejected at P2 rather than quietly degraded: skewed models with more than one slice, and motion without the averaging sliding band. The solver exits non-zero and says what is unsupported and how to get back to a working configuration. Nothing silently falls back to first order — a solver that downgrades the discretisation without telling you is worse than one that refuses.

And P1 results did not move: a nonlinear transient PMSM with motion, sliding band and circuit, and a 252-step induction machine with eddy currents and anti-periodic boundaries, both produce bit-identical solution files before and after the whole P2 feature. That is the regression standard worth holding a solver to.

Second-order elements are part of Nabla Core. Related: second-order time integration, far-field boundaries, and the validation dossier.