Fix Quadrature rule docs and syms

[termi-official]

Following the discussion in #1001 I just noticed that some references for the quadrature roles and the doc strings were not correct. This PR fixes both.

[KnutAM]

Nice that you found this!

Since we probably don't want to change the defaults in the future, are these reasonable defaults? Seems strange to switch at order = 8?

(Side-note: Also looks like dispatching on a type would be nicer as @fredrikekre suggested, but that would be a different PR)

[termi-official]

Nice that you found this!

Thanks, took me a while to realize too tbh.

Since we probably don't want to change the defaults in the future, are these reasonable defaults? Seems strange to switch at order = 8?

Actually I am really not sure about the rule or why it works in the in first place, because it comes with negative weights. I am pretty sure this can cause some trouble, but I cannot pin down why. See e.g. https://www.math.ntnu.no/emner/TMA4130/2021h/lectures/GaussQuadrature.pdf . But on the other hand rules with positive weights have more points in general.

(Side-note: Also looks like dispatching on a type would be nicer as @fredrikekre suggested, but that would be a different PR)

I am not super sure about this. I think it would also be nice to provide extensions which add new quadrature rules.

[codecov]

Codecov Report

Attention: Patch coverage is 93.50649% with 5 lines in your changes missing coverage. Please review.

Project coverage is 93.68%. Comparing base (2c8b751) to head (6b3ff91).
Report is 5 commits behind head on master.

Files	Patch %	Lines
src/Quadrature/gaussquad_tet_table.jl	89.47%	4 Missing ⚠️
src/Quadrature/gaussquad_tri_table.jl	87.50%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1007      +/-   ##
==========================================
+ Coverage   93.20%   93.68%   +0.47%     
==========================================
  Files          38       39       +1     
  Lines        5816     5936     +120     
==========================================
+ Hits         5421     5561     +140     
+ Misses        395      375      -20     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[fredrikekre]

The triangle ones are from #514 (comment) (unless you found it already).

[termi-official]

The triangle ones are from #514 (comment) (unless you found it already).

Yes, and going into the basics source code tracking down the call graph manually I arrived at a call to Gauss-Jacobi generators.

[termi-official]

Someone should verify this independently and then it should be good to go.

[KnutAM]

OK, started looking through this, and as I understand, we can skip the :jinyun type, as this is covered by :keast.

I've checked that :jinyn matches the reference so far (except for n=4 where it matches :keast), but don't have more time now.

[KnutAM]

Move to above _get_gaussjacobi_tridata and add actual ref to Gauss-Jacobi instead of the fenics link?

[termi-official]

I really struggle with finding the reference here. From what I understand it should be https://www.degruyter.com/document/doi/10.1515/crll.1826.1.301/html but the formula is slightly different ($w(x)=(x-1)^\alpha$?) from the formula I find in modern literature on ($w(x)=(x-1)^\alpha (x+1)^\beta$, e.g. https://epubs.siam.org/doi/pdf/10.1137/120889873).

[KnutAM]

But these are references for line intervals? Or do I miss something?

Seems like dunavant provides rules up to n=20. (https://onlinelibrary.wiley.com/doi/epdf/10.1002/nme.1620210612) (Although some are based on previous works, as noted on the bottom of page 1131).

Since higher than 7 were added in #514 and probably no one are saving quadrature data for such high orders, I think it could make sense to only use dunavant. But due to the different coordinate systems used in the papers and in the code, I find it hard to check where the points come from.
How did you check this to figure out that starting at n = 9 this is gauss_jacobi?

[termi-official]

But these are references for line intervals? Or do I miss something?

Yes and from my understanding you make a tensor product over the barycentric space to obtain the rule, but this is where I really got stuck.

Seems like dunavant provides rules up to n=20. (https://onlinelibrary.wiley.com/doi/epdf/10.1002/nme.1620210612) (Although some are based on previous works, as noted on the bottom of page 1131).

Since higher than 7 were added in #514 and probably no one are saving quadrature data for such high orders, I think it could make sense to only use dunavant. But due to the different coordinate systems used in the papers and in the code, I find it hard to check where the points come from.
How did you check this to figure out that starting at n = 9 this is gauss_jacobi?

Git blame and the form of quadrature points is very different (9,10,11 did not match the Dunavant paper and from the referenced PR I concluded that these must be the mentioned Gauss-Jacobi rules from FEniCS).

[KnutAM]

I think we should delete this option, since the current implementation is "keast" and we don't need to implement all possible quadrature rules. Of course, in the docs, we should acknowledge that keast orders 2 and 3 comes from jinyun.

[termi-official]

I would leave it for reference, because for me it was not that straight forward to figure out how to get the rules for order 1..3 from the Keast paper. This way it is directly obvious.

[KnutAM]

IMO, we should make this clear from the docs rather than adding more complexity/options to the code.

[KnutAM]

Checked keast 5 and keast positive 4 as well. They look ok (didn't check all points of course).

I think adding a test that checks that the defaults are actually used would be good.

[termi-official]

I think adding a test that checks that the defaults are actually used would be good.

Should we do this separately, as I think we have not yet closed the discussion on the defaults?

[KnutAM]

Should we do this separately, as I think we have not yet closed the discussion on the defaults?

Can be done in a separate PR, but no need to wait until the defaults discussion is closed I'm thinking that such a test is just to ensure that we don't change by mistake, and changing the test makes it clear that it was intended.