Redesign of quadrature rules
Currently there are some different ways where and how quadrature rules are defined and implemented. We have 1d rules, generated by Maxima scripts
xxx.mac from analytic multi-precision expressions, using a template for the
InitHelper class, provided in two ways, values as double and values as string. We have simplex quadrature rules, provided as double values and manually listed in the
SimplexQuadratureRule up to low order. And there are tensor-product rules constructed from lower dimensional rules, e.g. a tensor-product of 1d rules for all dimensions. These can be transformed to other geometries than cubes, resulting in very expensive rules (many points).
This is the current status. I propose a redesign of these rules: First of all, except for the 1d rules (and correspondingly the tensor-product rules based on these 1d rules) all rules must be provided as tabulated values of coordinates and quadrature weights. It is really difficult to write down analytic expressions for higher order rules on arbitrary geometry elements. Every few years a research group therefor publishes a new set of rules for a specific geometry, that may be better in the sense of symmetry of complexity, compared to older rules.
- To extend our quadrature rules frequently it would be nice to have a way to generate an implementation directly from tabulated values
- Providing high precision rules in text format allows to generate rules for any data type (up to the printed precision)
- There are tools available to extend the precision of quadrature rules (preserving symmetry constraints), based on tabulated values.
Thus, I propose to add a "library" of tabulated quadrature rules for all our geometry types and a converter (similar to the one implemented in the Maxima script for 1d) that generates rules in a uniform way. Additionally I propose to implement the rules using a string representation of the values and a string-to-float converter, like
std::stringstream to generate the concrete values on construction of the quadrature rule. Rules are implemented as singleton, thus a construction takes place only one and the additional performance consts are negligible.
The tabulated values should be given in a defined format, e.g.
# optional comment lines, to add reference to the publications where the rules were found dim order x0 y0 z0 weight0 x1 y1 z1 weight1
(see a set of rules in
utility/rules) 1d rules can be generated using the Maxima-Script as before, but saved in this tabulated format instead of generating c++ code directly.
I have already implemented a generator, based on Mako python templates (see
utility/make_rules.py, but it can simply be implemented in any language. The main difficulty is, that every publication provides values in a different format, for a different reference geometry, in barycentric coordinates or local coordinates, or .... Thatwhy I have additionally implemented a conversion script to read values given from a publication for one reference element and transform it to the Dune reference element and coordinate format.
Discussion from Dune Developer Meeting 2017
- all agree that the proposed changes are nice, as long as Simon does the changes
- the MR should be split up into smaller parts
- it is still up for dicussion, which intermediate format one should use
- it is a good idea to only use string representations and generate instances of float like data from this string representation
- in a future version it should be possible to honor further properties of a quadrature rule
- Cleanup of quadrature rules
- Return coefficients by cast from string representation
- Tabulation of quadrature rules
It was discussed whether to provide the tabulation in a way that other libraries can benefit from.
Suggestion: I create a new dune module
dune-quadrature that is responsible for the tabulation of quadrature rules and provides generator functions to produce C++ source files that can be included in