Hi there,

I have a finite element theory question (so likely for @tjhei or @bangerth, but anyone with FE experience is welcome to chime in):

By default we often use a Gauss quadrature with `FE_degree+1`

for our assemblers, because that can accurately integrate polynomials of order `2*n-1 = 2*FE_degree+1`

, and can therefore accurately represent terms like `phi_u * phi_u`

, which is of polynomial degree 2 * FE_degree. However we often use the same quadrature degrees in our postprocessors, where we simply want to integrate over the solution vector, which is of degree `FE_degree`

. If I am not mistake in order to integrate this accurately we would only need a Gauss quadrature of order `FE_degree/2 + 1`

. Is this correct? Because then we can reduce the quadrature order in many postprocessors (I am particularly looking at the depth average, because it is so expensive, and @maxrudolph needs it frequently in his models, every reduction in quadrature order would make it significantly cheaper).

Any advice is appreciated.

Rene

The way it is often taught is that a Gauss formula with q degrees is exact for

polynomials of degree 2q-1, and so to integrate finite element matrices like

the mass matrix (which has polynomial degree 2p), it is enough to use q=p+1.

In your case, you even have only polynomials of degree p, so q=roundup(p/2)

should be enough.

But this argument is not quite correct. It is only true for linear mappings

(on triangles, for example) because then the integrand really just remains a

polynomial of the same degree as on the reference cell. On the other hand, if

one is on quadrilaterals, or uses higher order mappings, then the actual

integrand one is integrating is no longer polynomial (but in reality some sort

of rational function). Consequently, Gauss quadrature of sufficiently high

order will no longer integrate *exactly*. Instead, what happens if you use

q=p+1 quadrature points is that the error one commits by replacing integration

by quadrature is *of the same order* as the finite element error, and so

sufficient to retain the overall convergence order. This is good enough. Using

a better quadrature formula would only reduce the constant in the error

(though likely not by very much in practice).

So your question is not so easy to answer. In practice, I suspect that if you

are integrating a quantity that is linear in the solution – say, the average

velocity/temperature/heat flux/… – then you are correct that a smaller

number of quadrature points will suffice. You won’t get the exact same answer,

but you shouldn’t get something that has a fundamentally different convergence

rate.

Best

W.

Thanks for the explanation!

I checked with the reduced quadrature order `q=roundup((p+1)/2)`

, and the difference in output is below 1e-3 for nearly all tests. There were some problems with exotic geometries (see https://github.com/geodynamics/aspect/pull/2910), but I suspect that is because of the geometries or model setup, not the postprocessing.

The depth averaging for `SphericalShell`

and `Box`

is nearly as accurate and significantly faster with the new version.

I checked with the reduced quadrature order |q=roundup((p+1)/2)|, and

the difference in output is below 1e-3 for nearly all tests.

Is this the *relative* difference? So is the difference in the third

significant digit of results? That really seems to be quite small, in

particular because almost all tests run on pretty coarse grids and the

error due to the discretization should still be quite large on these grids.

Yes that is relative difference on a coarse grid. Good enough I would say