Unify LagrangeSimplexFiniteElement across dimensions

The LagrangeSimplexFiniteElement contains 5 different code paths to implement the basis and its derivatives:

• A special implementation for order 0 in any dimension.
• A special implementation for order 1 in any dimension.
• An implementation for any order in dimension 1.
• An implementation for any order in dimension 2.
• An implementation for any order in dimension 3.

There are good reasons for the first two special cases. Unfortunately the ones for dim=1,2,3 are all implemented differently for no good reason (except history) makes the code hard to read and understand.

This patch unifies the latter three implementations. While there are still three branches, there is a very simple common structure making the now making the code much simpler to understand. The unified implementation basically follows the former evaluateFunction() implementation of the 3d case:

1. Rescale the simplex by k such that Lagrange points have integer coordinates.
2. Transform to barycentric coordinates in the scaled simplex.
3. Loop over Lagrange points in lexicographic order.
4. Compute basis functions as products of univariate Lagrange-polynomials wrt. barycentric coordinates.

This approach is now also used for evaluateJacobian() and first order derivatives in partial() extending the support of the latter to dim=1 and dim=2. arbitrary order partial derivatives in partial() for dim=1,2,3-

For now the special implementation of second order partial derivatives in 2d is kept, but could be updated accordingly.

This includes the cleaned up test from !293 (merged). Hence the latter should be merged before.

Additionally to the improved readability and added support for first order partial derivatives, the code is also shorter. While this also improves the performance (*), this is IMO of minor significance, because it doesn't matter when evaluations are cached.

(*) To check the performance I used this benchmark provided by @simon.praetorius which measures testFE(). It turns out that the performance is essentially identical in 3d but significantly improved in 1d and 2d. Notice that for 1d one needs to increase the number of test runs (e.g. by 100) to avoid the the differences are lost in the measurement perturbations. One may debate if this benchmark is significant, because it does much more work than done in actual computations. However, since only the bare computational code changed, one may expect that the performance improvement is even stronger if the test boilerplate is removed.

@oliver.sander As far as I can see, you're the one who changed the implementation last (but 6 years ago). Can you have a look?

Would this also solve the issue with the missing implementation of partial for 1d? see my comment 63a48974 (comment 144511)

As described above the implementation uses the barycentric decomposition L_{i}(x_0,...,x_{d-1}) = \prod_{j=0}^d \hat{L}_{i_j}(x_j) where x_d = 1-\sum_{j=0}^{d-1}x_j of the d-variate Lagrange polynomial L_{i} where the \hat{L}_{i_j} are univariate Lagrange polynomials. This can even be exploited further by precomputing all (d+1)(k+1) possible Lagrange polynomials \hat{L}_{i_j}(x_j).

In the current code this amounts int first storing the values lagrangeFactor(z[j], i) for z=0,...,d and i=0,...,k in a table and then simply computing all admissible products. Using this the benchmark performs even better, but I only implemented it for evaluateFunction() so far. But this should probably be tested for all orders indenpendently.

I think this MR is a very good idea. Simplifying the code shows better possible optimizations and might be less error prone.

A few improvements that seem necessary in order to support other range and domain types than double.

added 1 commit

• 8130945a - Precompute barycentric factors in LagrangeSimplexBasis

Compare with previous version

added 1 commit

• a598ce78 - Fix type in second order partial derivatives

Compare with previous version

resolved all threads

added 1 commit

• 218ab15d - Fix type in second order partial derivatives

Compare with previous version

I just pushed an updated version that

• precomputes the univariate Lagrange polynomials in barycentric coordinates,
• adds the formula for the basis functions to the doxygen documentation because this makes it easier to understand the code later on,
• includes the two return type fixes mentioned by @simon.praetorius in unrelated code branches.
• implements partial derivatives of arbitrary order (making the fixes to the old code paths obsolete).

This made the code again simpler and faster (in particular for higher order in 3d).

I now tested this by benchmarking testFE() independently for dim=1,2,3 and order 0 to 7. The new implementation is never slower and sometimes significantly faster (up to factor 6 and 3 for order 7 in 2d and 3d, respectivly). But most importantly, it's much better understandable.

added 7 commits

• 218ab15d...073573ea - 2 commits from branch master
• 68a798e9 - [test] Cleanup lagrangeshapefunctiontest
• 03a2b080 - Unify LagrangeSimplexFiniteElement across dimensions
• f1a9769a - Precompute barycentric factors in LagrangeSimplexBasis
• b1a6aa43 - Implement arbitrary order partial derivatives in LagrangeSimplexLFE
• e92dd53c - [test] Test 2st order partial() for LagrangeSimplexFiniteElement

Compare with previous version

changed the description

resolved all threads

Other types than double are still not tested, but this needs some adjustment to the testlocalfe.hh