Custom integrals fail for higher order derivatives

[Anders Logg (Chalmers)]

Code generated for the following example:

a = inner(div(grad(u)), div(grad(v)))*dc

fails to compile (on the C++ side) due to screwed up FE tables:

error: 
redefinition of 'FE0_C0_D11'
std::vector<std::vector<double> > FE0_C0_D11(num_quadrature_points);

[Prof Garth Wells]

Where is the std::vector<std::vector<double> >? In the generated code? std::vector<std::vector<double> > does come with some overhead because of the multiple memory allocations.

[Martin Sandve Alnæs]

I would expect that to dominate the cost if it's in tabulate_tensor.

[Anders Logg (Chalmers)]

It is used in place of the static arrays used for normal quadrature code, for example:

static const double FE1_C0_D10[6][15]

For custom_integral, the number of quadrature points is not known at compile time so it needs to be dynamics (or of a large fixed max size):

// Create table FE0_C0_D10 for basis function derivatives on all cells
std::vector<std::vector<double> > FE0_C0_D10(num_quadrature_points);
for (std::size_t ip = 0; ip < num_quadrature_points; ip++)
    FE0_C0_D10[ip].resize(15);

Yes, this is likely expensive and can be improved by flattening the arrays, but there are already something like 3 levels of offset computing (basis function index, component, derivative) so I didn't want to add a 4th until I knew things worked ok. Plus the code generation relies on reusing the standard code generation for quadrature as much as possible, hence a similar data structure.

[Anders Logg (Chalmers)]

Fix pushed to next

[Prof Garth Wells]

We could consider std::array from c++11, or boost::multi_array.