Making quadrature_rule a public method on relevant element groups

[thomasgibson]

This is based on a discussion in inducer/grudge#62.

Currently NodalElementGroupBase has a public user-facing method weights. This is used internally in grudge when integrating on these elements, using the unit_nodes as quadrature points.

Should we move _quadrature_rule (found in here and similar places) out and add quadrature_rule as a public method to NodalElementGroupBase. This would return a modepy.Quadrature object that can be used for reasoning about quadrature accuracy (and integration obviously 😄 ).

I am in favor of this. However, we have a problem that we need to address.

Problem:
As briefly described in this issue, some element groups do not have a concrete subclass of modepy.Quadrature (e.g. JacobiGaussQuadrature, LegendreGaussQuadrature, etc.) As a concrete example, consider subclasses of _MassMatrixQuadratureElementGroup. Grudge uses the weights method via the mass matrix and the interpolation nodes (unit_nodes) to compute using "mass matrix quadrature" (using @inducer's terminology here). It's easy enough to create a modepy.Quadrature object with the relevant weights/nodes, BUT when we have code that wants to check the order of accuracy via quad.exact_to, we run into problems.

Quick summary:
Creating a generic Quadrature object as a bag of nodes and weights lacks sufficient information to determine what exact_to needs to be. Currently, you can create a generic quadrature rule, but you cannot ask quad.exact_to (because it's NOT defined!)

If we want quadrature_rule to be a public-facing method, we should resolve the problem.

[inducer]

Thanks for writing up a summary of this! Suppose we go with the direction that inducer/modepy#33 seems to be leaning (?) (i.e. allow feeding exact_to to the Quadrature constructor), I claim that the resulting "mass matrix quadrature" is exact to the order of the reference polynomial space, by definition of the mass matrix as long as:

• the DOFs are nodal, and
• 1 is represented exactly, by a DOF vector of all ones.

[inducer]

(Opinions?) cc @alexfikl

[thomasgibson]

I claim that the resulting "mass matrix quadrature" is exact to the order of the reference polynomial space, by definition of the mass matrix as long as:

• the DOFs are nodal, and
• 1 is represented exactly, by a DOF vector of all ones.

I agree with this. This would allow quadrature_rule for the derived classes to return mp.Quadrature(self.weights, self.unit_nodes, exact_to=self.order), correct?

Though is this not well-conditioned, right? You mentioned at one point that this is equivalent to Newton-Cotes, which is... problematic for high order 😅

[inducer]

I agree with this. This would allow quadrature_rule for the derived classes to return mp.Quadrature(self.weights, self.unit_nodes, exact_to=self.order), correct?

Yep. At least as long as the polynomial space contains P^N (see the current, kinda hokey, definition of exact_to).

Though is this not well-conditioned, right? You mentioned at one point that this is equivalent to Newton-Cotes, which is... problematic for high order

Newton-Cotes with equispaced nodes is obviously terrible. Newton-Cotes with edge-clustered nodes is not that bad. (After all, Newton-Cotes with Legendre nodes is just Gaussian quadrature.)

Here's a demo: https://gist.github.com/inducer/19eb49ff7702d4bf5b12477208540731

[inducer]

(Note how even the weights stay positive in the demo.)

[inducer]

(But how the rule has obviously lost the 2*n accuracy of Gauss.)

[thomasgibson]

Here's a demo: https://gist.github.com/inducer/19eb49ff7702d4bf5b12477208540731

Very nice!

I've also gone ahead and started working on this issue: 44606a7