Calculation of surface terms could be faster

We could basically replace
https://github.com/trixi-framework/Trixi.jl/blob/b29bf00a7698d43160b393f739e1f8b5e37e057b/src/solvers/dgsem_tree/dg_3d.jl#L1137-L1165
by
```julia
@threaded for element in eachelement(dg, cache)
for m in eachnode(dg), l in eachnode(dg)
for v in eachvariable(equations)
# surface at -x
du[v, 1, l, m, element] = (
du[v, 1, l, m, element] - surface_flux_values[v, l, m, 1, element] * factor_1)

# surface at +x
du[v, nnodes(dg), l, m, element] = (
du[v, nnodes(dg), l, m, element] + surface_flux_values[v, l, m, 2, element] * factor_2)

# surface at -y
du[v, l, 1, m, element] = (
du[v, l, 1, m, element] - surface_flux_values[v, l, m, 3, element] * factor_1)

# surface at +y
du[v, l, nnodes(dg), m, element] = (
du[v, l, nnodes(dg), m, element] + surface_flux_values[v, l, m, 4, element] * factor_2)

# surface at -z
du[v, l, m, 1, element] = (
du[v, l, m, 1, element] - surface_flux_values[v, l, m, 5, element] * factor_1)

# surface at +z
du[v, l, m, nnodes(dg), element] = (
du[v, l, m, nnodes(dg), element] + surface_flux_values[v, l, m, 6, element] * factor_2)
end
end
end
```
to increase the performance for some common setups. For example, this change gives a 2x speedup of the surface calculation for
```julia
julia> trixi_include("examples/tree_3d_dgsem/elixir_euler_ec.jl",
volume_flux=flux_shima_etal_turbo, surface_flux=flux_shima_etal, tspan=(0.0, 10.0), polydeg=3)
```
This is measurable since the simple surface stuff takes more than 10% of the total runtime
```
──────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
────────────────────── ───────────────────────
Tot / % measured: 1.47s / 84.5% 16.9MiB / 92.3%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────
rhs! 1.19k 1.15s 93.1% 968μs 456KiB 2.85% 392B
volume integral 1.19k 481ms 38.8% 404μs 0.00B 0.00% 0.00B
interface flux 1.19k 297ms 24.0% 249μs 0.00B 0.00% 0.00B
surface integral 1.19k 153ms 12.3% 128μs 0.00B 0.00% 0.00B
prolong2interfaces 1.19k 149ms 12.1% 125μs 0.00B 0.00% 0.00B
reset ∂u/∂t 1.19k 33.7ms 2.72% 28.3μs 0.00B 0.00% 0.00B
Jacobian 1.19k 33.5ms 2.70% 28.1μs 0.00B 0.00% 0.00B
~rhs!~ 1.19k 4.79ms 0.39% 4.02μs 456KiB 2.85% 392B
prolong2boundaries 1.19k 256μs 0.02% 215ns 0.00B 0.00% 0.00B
prolong2mortars 1.19k 134μs 0.01% 112ns 0.00B 0.00% 0.00B
mortar flux 1.19k 107μs 0.01% 89.8ns 0.00B 0.00% 0.00B
boundary flux 1.19k 74.3μs 0.01% 62.4ns 0.00B 0.00% 0.00B
source terms 1.19k 42.9μs 0.00% 36.1ns 0.00B 0.00% 0.00B
calculate dt 239 37.1ms 2.99% 155μs 0.00B 0.00% 0.00B
analyze solution 4 35.7ms 2.88% 8.93ms 70.7KiB 0.44% 17.7KiB
I/O 5 13.2ms 1.07% 2.64ms 15.1MiB 96.7% 3.02MiB
save solution 4 12.2ms 0.99% 3.06ms 15.0MiB 96.4% 3.76MiB
~I/O~ 5 915μs 0.07% 183μs 41.9KiB 0.26% 8.38KiB
get element variables 4 39.7μs 0.00% 9.94μs 8.56KiB 0.05% 2.14KiB
save mesh 4 266ns 0.00% 66.5ns 0.00B 0.00% 0.00B
──────────────────────────────────────────────────────────────────────────────────
```

However, a change like this will need quite some discussion, see https://github.com/trixi-framework/Trixi.jl/pull/991#pullrequestreview-819957996