We can now specify an integration-rule for the SymbolicBFI, it works like this:
ir = IntegrationRule(points = [(0,0), (1,0), (0,1)], weights = [1/6, 1/6, 1/6] ) a += SymbolicBFI (u*v, intrule=ir)
A complete example performing mass lumping is in the attachment.