Add Vioreanu-Rokhlin quadratures #5
Conversation
… does not exactly integrate correct degree
|
I've added alternative quadratures for the problematic 9th and 11th order rules, taken from the `modepy' package and rescaled for our reference simplex. The quadratures are nearly identical and neither correctly integrate the polynomials as desired. Perhaps it is a numerical stability thing to do with our integrating on [0, 0], [0, 1], [1, 0] rather than over a domain symmetric about the origin? Not much more I can do here other than reduce the expected accuracies in the tests for those quadratures so tests succeed. (Maybe we should avoid using those orders, too). |
- Table 5.1 in Vioreanu-Rokhlin 2014 SISC is incorrect for three specific quadratures. - Also, test all degrees i + j <= p, for both VR and Gauss.
|
The author of this PR, tanderson92, is not an activated member of this organization on Codecov. |
Actual functionality is a TODO for reference_interpolation.jl.
test/reference_integration_test.jl
Outdated
| @testset "Vioreanu-Rokhlin quad on triangle" begin | ||
| d = Inti.ReferenceTriangle() | ||
| # exact value for x^a*y^b integrate over reference triangle | ||
| exa = (a, b) -> factorial(a) * factorial(b) / factorial(a + b + 2) | ||
| # check all quadrature implemented | ||
| orders = keys(Inti.TRIANGLE_VR_ORDER_TO_NPTS) | ||
| for p in orders | ||
| q = Inti.VioreanuRokhlin(; domain=d, order=p) | ||
| x, w = q() | ||
| @test Inti.domain(q) == d | ||
| @test all(qnode ∈ d for qnode in x) | ||
| for i in 0:p, j in 0:p | ||
| i + j > p && continue | ||
| @test Inti.integrate(x -> x[1]^i * x[2]^j, q) ≈ exa(i, j) | ||
| end | ||
| end | ||
| end |
There was a problem hiding this comment.
Maybe we should also check the quadrature does not integrate monomials of order
There was a problem hiding this comment.
I don't think that that makes sense; it isn't a property of being a quadrature of order p that monomials of order p+1 are not integrated. Indeed sometimes due to symmetries some of the monomials of order p+1 are integrated, but not all.
There was a problem hiding this comment.
If all monomials of order order function is ambiguous, since being of order
There was a problem hiding this comment.
Ok I have a commit that implements this. Funnily enough, Fejer is broken!
There was a problem hiding this comment.
Other than Fejer this is good to go from me.
There was a problem hiding this comment.
That is unfortunate. For some reason (symmetry?) Fejer is one order better for odd values of N. Should we fix this?
|
I just noticed this branch is failing on nighly. Usually not a big deal, but I wonder why... |
I noticed this too but the logs don't indicate that it has anything to do with our changes. I blame julia nightly. |
- FIXME: Fejer is broken!
tanderson92
changed the title
These quadratures are imported from a set of rules generated in accordance with the Vioreanu-Rokhlin 2014 SISC paper. The tables are by Gimbutas. They may or may not correspond to those in the "modepy" project.
TODO: