UFL like assembler framework

I played around with a simple UFL-like framework for defining assemblers. It turns out that this nicely interacts with the nested dune-functions bases.

To be more flexible while playing, it's in a separate project https://git.imp.fu-berlin.de/agnumpde/dune-fufem-forms so far. The whole implementation is essentially in the forms.hh header. Furthermore these are copies of the stokes-taylorhood, poisson-pq2, and poisson-mfem examples from dune-functions. For comparison they contain the new and old assembly code. In a nut shell you can do the following:

• Poisson with PQ2-discretization

auto v = testFunction(basis);
auto u = trialFunction(basis);
auto f = Coefficient(rhsFunction);

auto a = integrate(dot(grad(v), grad(u)));
auto b = integrate(f*v);

 // [...] Create global matrix, rhs, global dune-fufem assemblers
operatorAssembler.assembleBulk(matrixBackend, a);
functionalAssembler.assembleBulk(rhsBackend, b);

• Poisson problems with mixed RT/Lagrange-discretization

auto fluxBasis = Functions::subspaceBasis(basis, _0);
auto pressureBasis = Functions::subspaceBasis(basis, _1);

auto sigma = trialFunction(fluxBasis);
auto u = trialFunction(pressureBasis);
auto tau = testFunction(fluxBasis);
auto v = testFunction(pressureBasis);

auto A = integrate(dot(sigma, tau) + div(tau)*u + div(sigma)*v);

// [...]
operatorAssembler.assembleBulk(matrixBackend, A);

• Stokes problem with Taylor-Hood discretization

auto velocityBasis = Functions::subspaceBasis(taylorHoodBasis, _0);
auto pressureBasis = Functions::subspaceBasis(taylorHoodBasis, _1);

auto u = trialFunction(velocityBasis);
auto p = trialFunction(pressureBasis);
auto v = testFunction(velocityBasis);
auto q = testFunction(pressureBasis);

auto A = integrate( dot(grad(u), grad(v)) + div(u)*q + div(v)*p );
// [...]
operatorAssembler.assembleBulk(matrixBackend, A);

Non-constant coefficients should also work but are not tested so far. This is still WIP but works surprisingly nicely (with only about 700 loc). If there's interest, we can move this to dune-fufem. Current restrictions are:

• You need a patch in dune-common to make dot unambiguous.
• In a special case (dot-product on matrix-valued gradient for power of scalar space) this up to 40% slower, because it currently does not take advantage of the tensor-basis structure and multiplies many zeros (*). Writing the sum of scalar dot products explicitly avoids this.

Up to one exception (see (*) above) this is even faster than the hand-written implementations in the dune-functions examples.