Finite difference SBP methods

This is to track the progress of porting finite difference methods to Trixi. I propose to use my package [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl) to get the basic methods. Then, we could use the following approach (implemented in #617).
- [x] Pack an SBP derivative operator from SummationByPartsOperators.jl as `basis` in Trixi's `DG` struct
- [x] Organize methods written for general DG operators such that they can use these new basis types
- [x] Implement the weak-form volume integral in 2D for these SBP operators and the `TreeMesh` with conforming interfaces (no mortar, no AMR)

We shouldn't look at performance too much at this stage. It will probably be okay to just compute the fluxes and use `mul!`. However, that will require multiple `mul!`s for each component, `reinterpret`ing arrays to have `SVector`s as contents, or to use something like [StructArrays.jl](https://github.com/JuliaArrays/StructArrays.jl), see also #557.

Some limitations of this first approach and open questions are
- We will use the same number of nodes per coordinate direction
- Efficient parallelization will probably require new/more code. Will we have enough elements such that threading across elements is still optimal? Or will we have a few huge elements such that we also need threading inside of the elements?
- Truly multi-dimensional implementations can be more (cache-) efficient than dimension-by-dimension versions implemented at first
- How should we handle a single truly periodic block? The basic option above will give us one SBP element coupled via numerical fluxes. Do we also want to enable a single truly periodic FD block?
- Shall we also allow upwind SBP operators? If so, how can we do that efficiently?

This will be a basic proof-of-concept and allows us to see whether/how our current abstractions need to be changed. Afterwards, some possible next steps are
- [x] Rename volume/surface integrals to terms (#616)
- [ ] Implement purely periodic operators
- [x] Move code to appropriate places (cf. #542). We need more levels in our solver hierarchy.
- `DG`: Current sub-methods are FDSBP and DGSEM
- For each mesh type, we could have a `dg_common.jl`, a `dgsem_common.jl`, and a `fdsbp_common.jl`
- More specific methods can be implemented in `tree_2d_dgsem`, `tree_2d_fdsbp` etc.
- [ ] Add FDSBP methods to the benchmarks
- [ ] Add references
- [x] Add docs (including surface terms/integrals, e.g., https://github.com/trixi-framework/Trixi.jl/pull/617#discussion_r642628168)
- [x] Port stuff to 1D and 3D
- [ ] Extend `convergence_test` to different refinement types (number of elements, number of nodes per element)
- [ ] Implement split-form/flux differencing volume integrals. This could either be done in Trixi or moved (partially) to SummationByPartsOperators.jl, see #3187.
- [x] Implement curvilinear meshes #1773
- [ ] Adapt visualization to FDSBP methods
- [ ] Implement non-conforming meshes and AMR
- [ ] Implement shock-capturing stuff
- [ ] Create keyword constructors with default `mortar=nothing` for `DG` (https://github.com/trixi-framework/Trixi.jl/pull/617/files#r643664755) and add more explanation of `mortar` and the other terms (https://github.com/trixi-framework/Trixi.jl/pull/617/files#r644461598)
- [ ] Optimize the heck out of this stuff

Review the existing FDSBP and DG organization and the SummationByPartsOperators.jl integration. Choose one remaining item, such as purely periodic operators, benchmark coverage, visualization, or extending convergence_test, and first inspect the related entry point. Done means the selected item is implemented, integrated with the existing methods, and documented where applicable.

Written by the indexing model from the issue text.