ScalarEq: Implement modal Bspline1 DG evolution (!216)

Similar to !213

Same caveats apply, e.g. no HRSC working etc.

Continuation of !210

This adds a dispatch type `Bspline1` similar to `Bspline2`. Based on this we implement the necessary method to be able to use it in `DGElement`.

I also fixed a subtle bug in the quadrature rules I used in some of the `Bspline` methods.
That is, I accidentally used a LG rule with `Npts=2`, which corresponds to a `N=Npts-1=1` order rule that is exact for polynomials up to `2*N=2` order.
But for `Bspline2` we encounter fourth order polynomials, so we should be using `Npts=3`.

Also adds a new option to the `Mesh` parameters section:

* `kind = "modal_bspline1"`: Utilize a 1st order modal Bspline Ansatz for the DG evolution.

IntegrationTests: compactify test results summary (!214)

ScalarEq: Implement modal Bspline2 DG evolution (!213)

Similar to !207, but using the new 2nd order Bspline basis.

Currently supports only bare evolutions and no HRSC etc.

Also: This is a quasi-nodal Ansatz similar to the Bernstein Ansatz, i.e. it is not interpolating, but the modal coefficients are related to some sampling of the data on the grid.
This means that one has to sum the modal coefficients onto the grid to obtain the nodal values, which can be done with `interpolate` or a Vandermonde matrix.
Atm this is not consistently done in the ScalarEq project, hence, care needs to be examined when interpreting results!

This adds multiple things:
- A `Bspline` module for working with 2nd order Bsplines.
- Extends `DGElement` with the `kind=:modal_bspline2` option.
- Updates docstrings and refactors some tests.
- Adds a `LG` module for Legendre-Gauss quadrature.
- Adds a plotting script `plot/bspline.jl` to tinker with the `Bspline2` approximation.

Also adds a new option to the `Mesh` parameters section:
- `kind = "modal_bspline2"`: Utilize a 2nd order modal Bspline Ansatz for the DG evolution.

CI: try to reuse precompilation cache at least for IntegrationTests (!212)

dg1d: Fix convention for Bernstein polynomials (!211)

We now establish the convention that all basis functions are indexed by one.

Also had to update the `modal_bernstein_advection_sine` reftest data, because !201 skipped running CI on main and me cancelling some other CI runs.

We require that all PRs are rebased before merging, so running them twice does not help anything.
This should free up the CI queue.

ScalarEq: Implement modal Bernstein DG evolution (!207)

MWE to check consistency of operators.

This is not intended to work with any HRSC schemes or 2d etc.

dg1d: Rename `SpectralElement` to `DGElement` and generalize its purpose (!206)

Generalization should allow it represent DG operators that are not necessarily based on nodal approximations via Lagrange polynomials.

Atm we implement the `nodal_bernstein` version, which is actually pseudo-nodal and not yet tested. Will be done in a follow-up PR.

Along the way we also implement some operators for Bernstein and Legendre bases.

tests: add kw only_failed to `runtests` command (!205)

Unify interface for LGL and GLGL rules and update `SpectralElement` accordingly:
- most routine should just work with one argument `Npts`
- add more doc strings
- remove unused methods `linear_dependency_matrix` and `expansion_matrix`
- remove unused fields of `SpectralElement`

GRHD: fix another edge case where c2p_set_atmosphere_on_failure should overrule (!203)

dg1d: substep callbacks have to be evolved after updating the solution (!202)

dg1d: update `lgl` and `glgl` integration methods (!201)

- put the integrand as the first argument
- add integration methods that can work with a function argument
- add doc strings
- update usage of `integrate` functions in rest of code

GRHD: fix atmosphere freezing and evolution logic (!199)

Fixes the following issues with the current implementation:
- `update_atm_domain_of_dependence!` did not distinguish between points
  1. whose neighbors are also in the bulk,
  2. where one neighbor is on an interface,
  3. where the point itself is on an interface.

  These cases are already properly handled in 2d.
- `stop_atmosphere_evolution_*` methods did not account for the `c2p_freeze_atm` flag, which essentially prevented the atmosphere to ever evolve.
- `update_atm_domain_of_dependence!` was activated by setting `c2p_dynamical_atm`, but it should be controlled by `c2p_enforce_causal_atm`



---

Additional fixes with less impact:
- remove unused parameter `id_atmosphere_density` from `GRHD` section
- `c2p_set_atmosphere_on_failure` should really suppress any errors appearing in the cons2prim

I am surprised that I did not need this before.

This is intentionally a hacky implementation as I just need that for the thesis right now.

---

This also implements
- the total rest-mass computation for the GRHD project as a first application,
- adds helper functions for the `Output` struct to read all data at once.

Fixup to !193

GRHD: add tov symmetrization for all DG rhs' that use the spherical1d (!195)

dg1d: make `main` accept paths to output directories (!196)

dg1d: do not strip trailing newline chars as this might corrupt git patches (!194)