FEniCS/basix

Change quadrature rules

Closed this issue · 4 comments

jorgensd commented

In the package quadpy, there is a large collection of quadrature rules.
https://github.com/nschloe/quadpy
i would suggest replacing some of the existing rules with rules from here. The rules are stored as json files.

mscroggs commented

Progress towards this was made in #258

garth-wells commented

@jorgensd Which rules in particular?

jorgensd commented

Whenever we fall back to an mth order Gauss-Jacobi rule (https://github.com/FEniCS/basix/blob/main/cpp/basix/quadrature.cpp#L1736-L1759), we use:

  • Interval: (m+2) / 2
  • Quadrilateral ((m+2) / 2)**2
  • Hexahedron ((m+2) / 2)**3
  • Triangle: ((m+2) / 2)**2
  • Tetrahedron: ((m+2) / 2)**3

Meaning that for a 9th order quadrature rule on a tetrahedron, we have 125 quadrature points while Xiao-Gimbutas only uses 57 points (as far as I can tell). https://github.com/nschloe/quadpy/tree/main/src/quadpy/t3/_xiao_gimbutas
The goal would be to get all the default results to a relatively high order to use other points than Gauss-Jacobi.

Do you have anything to add @mscroggs ?

mscroggs commented

I think we have good schemes on intervals, quads, hexes, and triangles. Once we add some rules for tets up to a reasonable order, I think this is done