Generate shape functions

[Psirus]

Currently, we hardcode the shape functions and their derivatives, ending up with code that looks like this:

    //    shapeFunctions[13]=  0.25 *  plus_r * minus_s *   mid_t;
    derivativeShapeFunctions(13, 0) = 0.25 * minus_s * mid_t;
    derivativeShapeFunctions(13, 1) = -0.25 * plus_r * mid_t;
    derivativeShapeFunctions(13, 2) = -0.5 * t * plus_r * minus_s;

    //    shapeFunctions[14]=  0.25 *  plus_r *  plus_s *   mid_t;
    derivativeShapeFunctions(14, 0) = 0.25 * plus_s * mid_t;
    derivativeShapeFunctions(14, 1) = 0.25 * plus_r * mid_t;
    derivativeShapeFunctions(14, 2) = -0.5 * t * plus_r * plus_s;

    //    shapeFunctions[15]=  0.25 * minus_r *  plus_s *   mid_t;
    derivativeShapeFunctions(15, 0) = -0.25 * plus_s * mid_t;
    derivativeShapeFunctions(15, 1) = 0.25 * minus_r * mid_t;
    derivativeShapeFunctions(15, 2) = -0.5 * t * minus_r * plus_s;

I think it would be good to generate these functions automatically. This would lead to less code, and would probably be easier to extend to arbitrary order. The main problem is that we probably would want this to be done at compile-time. So any constexpr-wizards are welcome to chime in :grinning:

[Psirus]

So I wanted to try this, and I have actually managed to generate "arbitrary" order (for order 71, it takes about 7 seconds to compile^^) Lagrange polynomials at compile time. The code is here and is still quite a mess of failed tries and intermediate steps that I needed to get to grips with the template meta programming stuff. I left it in, because it might be helpful to understand what is going on.

a) It uses boost fusion, although maybe boost mpl would suffice (or even be better?) b) still inconsistent in terms of where the ::type goes, in the "creator" or the "caller" (not sure what the correct terminology is).

Not sure if I can generate the lines

shapeFunctions[0] = LagrangePolynomial<3, 0>::F(x);
shapeFunctions[1] = LagrangePolynomial<3, 1>::F(x);
shapeFunctions[2] = LagrangePolynomial<3, 2>::F(x);
shapeFunctions[3] = LagrangePolynomial<3, 3>::F(x);

at compile time. Suggestions welcome.

Edit: A much simplified version can now be found here. No more boost :)

[joergfunger]

The new version looks much simpler than the original boost version. I very much like the idea, since it would generalize to any dimension. The only remaining ingredient to make this work in general is the determination of the nodal positions for either the spectral elements or other higher order schemes (e.g. hierarchical basis functions). Up til now, this is not done at compile-time, but once at runtime.