Fast(er) paths for function evaluation in reference space

In https://github.com/Ferrite-FEM/FerriteViz.jl/pull/167 we evaluate very often the finite element solution with known reference space points, so its really part of the hot loop and I tried to squeeze out as much performance as possible with the robot and it is possible to optimize `FunctionValues` path for different cases

```
┌─────────────────────┬─────────────┬────────────────┬─────────────────────────────┬──────────┐
│ Case │ PointValues │ FunctionValues │ direct sum, memoized vector │ monomial │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q1 quad scalar │ 10.0 │ 6.7 │ 4.5 │ 7.5 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q1 quad vector (^2) │ 15.5 │ 11.5 │ 6.5 │ 7.4 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q2 quad scalar │ 27.1 │ 24.0 │ 27.0 │ 13.4 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q2 quad vector (^2) │ 55.2 │ 48.1 │ 25.3 │ 13.3 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ P1 tri scalar │ 9.9 │ 6.3 │ 2.7 │ 7.0 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ P2 tri scalar │ 12.4 │ 8.9 │ 5.6 │ 12.2 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ P2 tri vector (^2) │ 58.6 │ 55.9 │ 9.4 │ 12.4 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q1 hex scalar │ 40.2 │ 22.3 │ 32.4 │ 12.4 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q1 hex vector (^3) │ 128.8 │ 81.2 │ 33.0 │ 13.6 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q2 hex scalar │ 83.5 │ 69.4 │ 89.0 │ 52.3 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ Q2 hex vector (^3) │ 364.7 │ 237.7 │ 88.3 │ 56.5 │
├─────────────────────┼─────────────┼────────────────┼─────────────────────────────┼──────────┤
│ P2 tet vector (^3) │ 140.5 │ 97.5 │ 36.9 │ 32.6 │
└─────────────────────┴─────────────┴────────────────┴─────────────────────────────┴──────────┘

```

- we could memoize the shape function value in the vectorized case which saves in 3D quite some time, since the values are the same just for different components (memoized vector column)
- for certain interpolations we can write it as
$$u_h(\boldsymbol{\xi}) = \sum_{i} N_i(\boldsymbol{\xi}), u_i = \sum_{m} c_m, \boldsymbol{\xi}^{\boldsymbol{\alpha}_m}, \qquad \mathbf{c} = V^{-1} \mathbf{u}e, \quad V{jm} = \hat{\boldsymbol{\xi}}_j^{\boldsymbol{\alpha}_m}$$

with $\hat{\boldsymbol{\xi}}_j$ the interpolation's reference nodes and $\boldsymbol{\alpha}_m$ the monomial exponents. $V^{-1}$ is computed once per interpolation, $\mathbf{c}$ once per cell and solution, and each point costs one multiply-add per monomial.

Not sure if there are other applications that could benefit from this but the vectorized memoization stuff seems worth imo