Default quadrature rule for Tets contains negative weight

[termi-official]

https://github.com/Ferrite-FEM/Ferrite.jl/blob/98e2b674a8c6d740a779c022de31f73051cb353c/src/Quadrature/gaussquad_tet_table.jl#L29

Negative quadrature weights can cause all sorts of issues when integrating functions, which are outside the basis space they are designed for. While such rules can be more efficient on many problems, I do not think that we should provide a problematic rule as the default rule to users, such that users must explicitly make the decision to use a rule with negative weights if they want the additional performance. After a bit googling I think https://portwooddigital.com/2020/07/03/integration-points-with-negative-weight/ explains the issue, which can be really hard to track down.

Hence I would propose that we should provide, as the default, only rules with positive weights to not surprise new users. Furthermore I think it makes sense to spend a few lines on this in the docs.

[KnutAM]

Negative quadrature weights can cause all sorts of issues when integrating functions, which are outside the basis space they are designed for.

Do you here mean the basis space that the quadrature rule is designed for?

While I see that the why is hard to explain. Do you know of any references that report that his is an issue for the tet cases where we have negative weights currently?

The blogpost seems to be more focused on structural elements, and in particular says

These approaches lead to varying degrees of strange results, but are generally harmless because the negative integration weights are not at the ends of the element. However, there could be unexpected results with member loads.

I interpret this as saying that negative weights are only an issue at the end points.

The example deal with plastic hinges and shows what happens when such a hinge is outside the element, leading to a negative integration weight. But I don't understand how the fact that an input error in defining the plastic hinge which gives a negative weight relates to the properties of having negative weights in general.

[termi-official]

Negative quadrature weights can cause all sorts of issues when integrating functions, which are outside the basis space they are designed for.

Do you here mean the basis space that the quadrature rule is designed for?

Yes. E.g. for some of the rules the space is polynomials of order p on the reference geometry.

While I see that the why is hard to explain. Do you know of any references that report that his is an issue for the tet cases where we have negative weights currently?

No reference here at hand, but you can just take any problem with localized nonlinearity and positivity. For example you can center the nonlinearity (e.g. something like exp(-100(x-c)^2)) on the point with the negative weight and obtain entries with the wrong sign in your stiffness matrix.

[KnutAM]

Valid argument!

[termi-official]

I know it is a bit made up, and most of the users won't ever hit such a case. However, I am just worried that, if users hit such a case and are not aware of this issue, then I think it will be insanely difficult for them to track the problem down, because the vast majority of the users come with the expectation that the quadrature rule does the "correct integration".

[KnutAM]

Yes, that was my thought too but sufficient to warrant playing it safe (a ref to motivate in the docs would be nice of course, but apart from that - after trying to get a bit into this, I think moving forward with new defaults should be ok).

Without knowing, I would guess that we sacrifice some accuracy for standard cases by changing to positive order here if measured in the highest exact polynomial order per quadrature point. On the other hand, for non-polynomial functions I'm not sure. But surely there must be results for this in the literature for e.g. rational, exponential, and other functions?

Since people are using Ferrite in so many different ways, favoring safe over fast/accurate default seems reasonable to me in this case. I guess the only place this will commonly affect users is L2 projection on tetrahedrons.

[termi-official]

Without knowing, I would guess that we sacrifice some accuracy for standard cases by changing to positive order here if measured in the highest exact polynomial order per quadrature point. On the other hand, for non-polynomial functions I'm not sure.

I think regarding accurracy there should not be a significant difference. Also note that most of the quadrature rules out there are not designed to integrate full polynomial systems, but slight variations (e.g. Laguerre for \int p(x)exp(-x)dx ,...).

But surely there must be results for this in the literature for e.g. rational, exponential, and other functions?

I was also wondering if we have some rules which can integrate the function found in for e.g. nonlinear geometry exactly.

[fredrikekre]

I think it is fine to changes this later without call it breaking.

[KnutAM]

Reference to motivation why this is bad (to not loose in the Slack-hole): https://www.sciencedirect.com/science/article/pii/S0898122121002820?via%3Dihub