Custom quadrature rules, runtime quadrature

[michalhabera]

There are few code paths for runtime quadrature. Unfortunately, this is identified with a custom UFL measure, while I think it should be just stored in measure's metadata.

One example interface for custom and runtime quadrature could be

v * dx(metadata={"quadrature_degree": 2, "quadrature_scheme": "whatever_preimplemented"})
v * dx(metadata={"quadrature_scheme": "custom", "quadrature_points": [[1.0, 0.5], ...], "quadrature_weights": [1.0, ...]})
v * dx(metadata={"quadrature_scheme": "runtime"})

FFCx signature for runtime quadrature (custom integral type) needs updating.

[bleyerj]

I don't know if this is the correct place to discuss this but many applications could use quadrature rules using vertices as quadrature points e.g. the already implemented "vertex" scheme on triangles, with quadrature points matching the vertices of a P1 Lagrange element. Other examples, may include Gauss-Lobatto points on lines or their tensorial version for quads/hex up to any order and quadrature scheme with points matching the vertices of Pk Lagrange element. In order to avoid having to define explicitly the points using a "custom" scheme for this situation, one can imagine a predefined scheme such as: dx(metadata={"quadrature_degree": 2, "quadrature_scheme": "nodal"}) or something else with a better name. This would allow to define very easily lumped mass matrix since u * v* dx(metadata={"quadrature_degree": 2, "quadrature_scheme": "nodal"}) would automatically be diagonal for Lagrange elements.

[chrisrichardson]

Can we use the same signature as the other integrals, and pack the points/weights as coefficients?

[michalhabera]

@bleyerj Thanks for the comment. Yes, the right place. So any quadrature scheme "default"/"vertex"/... etc. which is implemented in FIAT for the specific cell should work out of the box. The work to implement new quadrature scheme into FIAT should be perpendicular. I agree it would be nice to have more pre-implemented in stock. I have started working in this concept, will make PRs very soon.

@chrisrichardson I don't know if packing into coefficients array is best. I think we should add one scalar_type * for once and all. There are many cases where this would be useful. In general, we'd need a const pointer version as well for performance. And the common signature is a must, there are cases where one generated tabulate_tensor routines will have some quadrature done runtime and some compile-time. This was not possible with the old code, but I am actively working on this.

[chrisrichardson]

I tried the custom integral with cutfem last year. I recall there being issues with basis derivative evaluation.

[augustjohansson]

Thank you Michal for bringing up this topic! For what is worth, I personally think the format

v * dx(metadata={"quadrature_scheme": "custom", "quadrature_points": [[1.0, 0.5], ...], "quadrature_weights": [1.0, ...]})

is intuitive when working with cutfem-style methods. Of course, one only wants these on some parts of the domain, so typically one would like to have

v * (dx(list of elements with standard quadrature) + dx(list of elements with custom quadrature, ...))

and something similar for ds. I'm sure you know all this, but maybe it'll be helpful for some.

I can't comment on appropriate signatures more than to say that for these methods that we call multimesh methods, we take in pairs of elements (as in DG) and custom quadrature.

[michalhabera]

@bleyerj See #261

[sclaus2]

@michalhabera @mscroggs @chrisrichardson

Hi,

I would like to revive this discussion again because I am very interested in porting the CutFEM library to FEniCS-X in the coming months. In CutFEM and related methods, quadrature rules can vary from element to element and even within one element (multiple sub-triangles per element for integration). This means that hundreds or thousands of different quadrature rules are needed and therefore I think that a tabulate_tensor function does need to be generated for custom integral measures that can take in a quadrature rule at run-time in the assembly loop. A rough outline of the generated code would be

void tabulate_tensor( double* restrict A, const double* restrict w,
                 const double* restrict c, const double* restrict coordinate_dofs,
                 const int* restrict entity_local_index,
                 const uint8_t* restrict quadrature_permutation,
                 const int* restrict num_quadrature_points,
                 const double* restrict quadrature_points,
                 const double* restrict quadrature_weights)
{
    //Instead of pre-computing the values for shape functions and their derivatives 
   // they need to be computed here at the given quadrature points

    ...
    // generate code here to obtain information on the element and number of derivatives
    ...
    // Tabulate values of shape functions and their derivatives at given points
    //TODO: expose basix functionality to C and include appropriate header files in generated files
    basix_element_tabulate(deriv_order, quadrature_points,tab_data)

    // FE* dimensions: [permutation][entities][points][dofs]
    static double FE4_C0_D10[perm][entities][num_quadrature_points][ndofs];

    //Copy over info from tab_data to FE*..

    //Quadrature loop
    for (int iq = 0; iq < num_quadrature_points; ++iq)
    {
        ...
    }
}

In summary, the main extensions I see that are needed for this tabulate_tensor function for custom measures are to expose some of the basix functionality to tabulate basis values and basis derivative values to the generated C-code and to replace the pre-computation step in the tabulate_tensor code generation with the generation of code for the computation of the basis values with basix function calls as outlined above.