Skip to content

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.

This is a follow up issue of issue #6 . See the merge request !27 (closed).

Update

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 discussion, 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

Workflow

  • Cleanup of quadrature rules (see !138 (merged))
  • Return coefficients by cast from string representation (see !175 (merged))
  • 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 dune-geometry directly.

Edited by Simon Praetorius