Improve plotting for 2D/3D simulations with `polydeg = 0` on the `TreeMesh` (finite volume case)

Basic support for polydeg = 0 on the TreeMesh was implemented in https://github.com/trixi-framework/Trixi.jl/pull/1489, including a specialization of the plotting routines using Plots.jl. It would be nice to extend this plotting support to higher dimensions.

For example, I currently get

julia> using Trixi, Plots

julia> equations = LinearScalarAdvectionEquation2D(1.0, 1.0)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation2D                                                                  │
│ ═══════════════════════════════                                                                  │
│ #variables: ……………………………………………………………… 1                                                           │
│ │ variable 1: ………………………………………………………… scalar                                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> initial_condition(x, t, equations) = sinpi(x[1]) * sinpi(x[2])
initial_condition (generic function with 1 method)

julia> solver = DGSEM(polydeg = 0, surface_flux = flux_lax_friedrichs)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64}                                                                                      │
│ ═══════════                                                                                      │
│ basis: …………………………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=0)                    │
│ mortar: ………………………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=0)                 │
│ surface integral: ……………………………………………… SurfaceIntegralWeakForm                                     │
│ │ surface flux: …………………………………………………… FluxLaxFriedrichs(max_abs_speed)                            │
│ volume integral: ………………………………………………… VolumeIntegralWeakForm                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), n_cells_max = 10^5, initial_refinement_level = 2, periodicity = true)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2, Float64}}                                                        │
│ ═════════════════════════════════════════                                                        │
│ center: ………………………………………………………………………… [0.0, 0.0]                                                  │
│ length: ………………………………………………………………………… 2.0                                                         │
│ periodicity: …………………………………………………………… (true, true)                                                │
│ current #cells: …………………………………………………… 21                                                          │
│ #leaf-cells: …………………………………………………………… 16                                                          │
│ maximum #cells: …………………………………………………… 100000                                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; boundary_conditions = boundary_condition_periodic)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic                                                                     │
│ ════════════════════════════                                                                     │
│ #spatial dimensions: ……………………………………… 2                                                           │
│ mesh: ……………………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2, Float64}} with length 21    │
│ equations: ………………………………………………………………… LinearScalarAdvectionEquation2D                             │
│ initial condition: …………………………………………… initial_condition                                           │
│ boundary conditions: ……………………………………… Trixi.BoundaryConditionPeriodic                             │
│ source terms: ………………………………………………………… nothing                                                     │
│ solver: ………………………………………………………………………… DG                                                          │
│ total #DOFs per field: ………………………………… 16                                                          │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> ode = semidiscretize(semi, (0.0, 10.0));

julia> ode.u0
16-element Vector{Float64}:
  0.4999999999999999
  0.4999999999999999
  0.4999999999999999
  0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
 -0.4999999999999999
  0.4999999999999999
  0.4999999999999999
  0.4999999999999999
  0.4999999999999999

julia> plot(ode.u0, semi)

We clearly see that some interpolation has been performed although the numerical solution should have only two values.