Gauss quadrature in finite elements: how many points

An element stiffness matrix is an integral over the element:

Nobody evaluates it analytically. It is computed by sampling the integrand at a few points and taking a weighted sum, and the choice of how many points is one of the few decisions in a finite element program that changes the answer qualitatively rather than by a percentage.

Every number in this article is recomputed from its inputs by a script that runs on each build of this site.

The rule and its one propertyPermalink to “The rule and its one property”#

Gauss-Legendre quadrature on $[-1, 1]$ approximates

with abscissae and weights chosen so that the rule is exact for polynomials of the highest degree possible. Both the $n$ points and the $n$ weights are free, so $2n$ conditions can be satisfied, and the rule is exact to degree $2n - 1$.

That single property — exact to $2n-1$ — is the whole design rule. Everything below follows from it.

What under-integration costs on a plain integralPermalink to “What under-integration costs on a plain integral”#

Take $\int_{-1}^{1}\xi^4\,\mathrm{d}\xi$, whose exact value is $2/5 = 0.4$:

The one-point rule returns exactly zero, because the integrand vanishes at $\xi = 0$ and the rule sees nothing else. This is not a rounding matter. The rule is blind to a function that happens to be zero where it looks, and that blindness is the mechanism behind everything in the next section.

Three points integrate the quartic exactly because $4 \le 2(3) - 1 = 5$. Move to $\xi^6$ and three points give 0.240000 against $2/7$; four are needed, since $6 \le 2(4)-1 = 7$.

Spurious modes: the same blindness, in an elementPermalink to “Spurious modes: the same blindness, in an element”#

A plane element with $m$ nodes has $2m$ degrees of freedom, of which three are rigid body motions. The remaining $2m - 3$ are deformation modes, and each must be resisted by some stiffness or the element is singular in that direction. Each integration point supplies three independent constraints in plane analysis, one for each in-plane stress component.

Counting the two against each other:

A spurious mode is a displacement pattern that produces zero strain at every sampling point while producing real strain everywhere else. The element stores no energy in it, so the assembled stiffness matrix has a zero eigenvalue, and the mode is free to grow without bound. In the four-node quadrilateral these are the two hourglass modes, so called from the shape a row of them makes.

The counting is necessary but not sufficient — it says how many modes can be spurious, not that they are — and in a real mesh boundary conditions often suppress them. “Often” is the problem. A mode that is restrained in one mesh and free in the next, with no warning either way, is precisely the kind of behaviour a solver should not have.

When exact integration is impossiblePermalink to “When exact integration is impossible”#

For an undistorted rectangular element the integrand is a polynomial and there is a rule that integrates it exactly. For a distorted element there is not, and no number of points will produce one.

The reason is the Jacobian. Mapping from the reference square to a general quadrilateral introduces a factor $\det\mathbf{J}$ that varies over the element, and the strain-displacement matrix carries $\mathbf{J}^{-1}$. The integrand therefore contains a rational function — polynomials divided by $\det\mathbf{J}$ — and Gauss-Legendre is a polynomial rule.

So on any mesh with distorted elements, “full integration” is a convention, not a statement of exactness: it means the rule that would be exact if the element were a parallelogram. This is worth knowing because it is one of several reasons why two programs using nominally the same element can disagree in the fourth digit, and why a benchmark reproduced to 0.03% needs its mesh stated along with its answer.

Why anybody under-integrates on purposePermalink to “Why anybody under-integrates on purpose”#

Given the above, reduced integration looks indefensible. It is used anyway, for two reasons.

The first is cost. The number of points goes as $n^2$ in two dimensions and $n^3$ in three, and the constitutive routine is called at every point of every element at every iteration. Dropping a 3×3 rule to 2×2 removes five of every nine constitutive evaluations from a nonlinear analysis.

The second is better: full integration makes elements too stiff in near-incompressible conditions. The volumetric constraint is enforced at every integration point, so a fully integrated element carries more constraints than it has freedoms to satisfy them with, and it locks — returning collapse loads that are systematically high. Reducing the number of points reduces the number of constraints, which is exactly the cure, and exactly the disease described three sections above.

The honest resolutions are selective integration, which under-integrates only the volumetric term while integrating the deviatoric term fully; stabilisation, which adds an artificial stiffness targeted at the hourglass modes and nothing else; or higher-order elements, which have enough freedoms that full integration does not lock them. The last is why the fifteen-node triangle exists and why element choice is a modelling decision rather than a preference: with a lower-order element the computed limit load is wrong in a known direction, and no mesh refinement removes it.

That last point is worth stating plainly, because it is the practical consequence of all of the above: element type, integration rule and mesh are not three independent settings. They are one decision, and a result quoted without all three is not reproducible — which is the same argument that a verification suite rests on.