Support for arbitrary quadrature rules?

[Andrew Taber]

I'm wondering if FEniCS supports arbitrary quadrature rules (i.e. user-provided points and weights)? This issue (https://bitbucket.org/fenics-project/dolfin/issues/33/allow-for-custom-quadrature-rules) seems to suggest that custom quadrature rules are incorporated, but it looks like it's the degree of the quadrature rule that's defined at runtime, not the points/weights themselves.

Thanks!

[Jan Blechta]

There are three quadrature schemes supported: "default", "canonical", "vertex", see https://bitbucket.org/fenics-project/fiat/src/faffacc2cec53e1500c50272cb6594b6d06a51f1/FIAT/quadrature_schemes.py?at=master&fileviewer=file-view-default.

Adding support for arbitrary quadrature schemes is certainly possible and requires some work in FIAT, UFL, FFC and DOLFIN. Interesting thing is that the most of the work can be done in pure Python packages: FIAT, UFL, FFC, so this is low-hanging fruit but the approach needs to be discussed first on #development channel on FEniCS Slack.

[Andrew Taber]

Do you think using QuadratureElement or manually instantiating a different QuadratureRule when I'm building the forms is a fruitful avenue? It certainly looks like QuadratureElement was supposed to afford arbitrary quadrature functionality. Thanks!

[Jan Blechta]

No. Purpose of quadrature element is different. It is described in one chapter of the FEniCS book.

Quadrature element will only be able to use predefined quadrature schemes. It does not bring custom schemes.

[Lawrence Mitchell]

FWIW, with the tsfc representation, you can pass a QuadratureRule object: see for example https://github.com/firedrakeproject/tsfc/blob/master/tests/test_underintegration.py#L39. Although it's possible this does bad things for UFL's metadata handling.

[Jan Blechta]

Yes, probabably won't pass through DOLFIN/FFC parameter pipeline, but could work through integral metadata.

[Jan Blechta]

This is how one can monkey patch FIAT and FFC to do arbitrary quadrature. This is not intended for deploying by users but rather for demonstrating that this really a lot of music for just a pittance, or a low-hanging fruit for new developers.

Be careful with the cache! Recommend rm -rf tmp; DIJITSO_CACHE_DIR=$PWD/tmp python3 .... HIC SVNT LEONES!

from numpy import array

from FIAT.reference_element import UFCTriangle, UFCTetrahedron
from FIAT.quadrature import QuadratureRule
from FIAT.quadrature_schemes import create_quadrature
from ffc.analysis import _autoselect_quadrature_rule

def create_quadrature_monkey_patched(ref_el, degree, scheme="default"):
    """Monkey patched FIAT.quadrature_schemes.create_quadrature()
    for implementing custom scheme"""
    # Our "special" scheme
    if scheme == "bogus":
        if isinstance(ref_el, UFCTriangle) and degree == 42:
            x = array([[0.0, 0.0],
                       [1.0, 0.0],
                       [0.0, 1.0]])
            w = array([100.0, 200.0, -300.0])
            return QuadratureRule(ref_el, x, w)
        raise NotImplementedError("Scheme {} of degree {} on {} not implemented"
            .format(scheme, degree, ref_el))

    # Fallback to FIAT's normal operation
    return create_quadrature(ref_el, degree, scheme=scheme)

def _autoselect_quadrature_rule_monkey_patched(*args, **kwargs):
    """Monkey patched ffc.analysis._autoselect_quadrature_rule()
    preventing FFC to complain about non-existing quad scheme"""
    try:
        return _autoselect_quadrature_rule(*args, **kwargs)
    except Exception:
        integral_metadata = args[0]
        qr = integral_metadata["quadrature_rule"]
        return qr

# Monkey patch FIAT quadrature scheme generator
import FIAT
FIAT.create_quadrature = create_quadrature_monkey_patched

# Monkey patch FFC scheme autoselection mechanism
import ffc.analysis
ffc.analysis._autoselect_quadrature_rule = _autoselect_quadrature_rule_monkey_patched


# Standard FEniCS program
from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 2)
class X(Expression):
    def eval(self, values, x):
        values[0] = x[0]
bc = DirichletBC(V, X(degree=1), lambda x, b: b)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx(scheme="bogus", degree=42)  # Use custom quadrature
L = Constant(0)*v*dx
u = Function(V)
solve(a == L, u, bc)
import matplotlib.pyplot as plt
plot(u)
plt.show()

[Andrew Taber]

Following up on this issue after the discussion in the #development channel.

I should have been more specific in the body of the issue. My application is unfitted FEM, and therefore I need runtime per-element custom quadrature rules (each unfitted element gets its own quadrature rule in order to accurately integrate the intersection between the cell and the geometry).

I'm willing to put significant development time into this to make it happen. Ideally, whatever interface which solves my problem could also be used by other projects, such as as the CutFEM project.

What about an interface that somehow makes per-element queries to a dictionary of quadrature rules (either in-memory or on-disk, due to the potentially large amount of data) during assembly? This interface is agnostic to the method by which quadrature rules are generated (in my case I'm doing moment-fitting, CutFEM and MultiMesh have their own procedures to generate quadrature rules, etc).

[Nico Schlömer]

[slightly off-topic] Perhaps the existence of quadpy provides more motivation for this to be implemented.

[Jan Blechta]

@nschloe nice!

[Jan Blechta]

Closing this for inactivity but mainly because the issue in question is not in any way a DOLFIN issues. The issue concerns FIAT and UFL. (It turns out that UFL could carry FIAT (or quadpy?) cargo (which would be hashable) to describe quadrature - without making UFL depend on FIAT.)

[Jan Blechta]

New developers interested in this, please get in touch on FEniCS Slack (#development) channel.