diff --git a/src/discretization.jl b/src/discretization.jl
new file mode 100644
index 0000000000000000000000000000000000000000..bb81ed76d280a396cf4f68e42127233a91d5920b
--- /dev/null
+++ b/src/discretization.jl
@@ -0,0 +1,146 @@
+Base.@kwdef struct Discretization{
+ T_Grid<:AbstractGrid,
+ T_Cache<:AbstractCache,
+ T_Equation<:AbstractEquation,
+ T_HRSC<:Maybe{AbstractHRSC}}
+
+ grid::T_Grid
+ cache::T_Cache
+ equation::T_Equation
+ hrsc::T_HRSC = nothing
+end
+
+
+struct ParameterDict <: AbstractDict
+ provided::FlatDict{String,Any}
+ consumed::FlatDict{String,Any}
+ defaulted::FlatDict{String,Any}
+
+ function ParameterDict(d::AbstractDict)
+ flat_d = FlatDict(d)
+ consumed = FlatDict()
+ return new(flat_d, consumed)
+ end
+end
+
+
+function Base.getindex(pdb::ParameterDict, key)
+ val = get(pdb.provided, key, nothing)
+ if !isnothing(val)
+ if !haskey(pdb.consumed, val)
+ pdb.consumed[key] = val
+ end
+ return val
+ end
+ throw()
+end
+
+
+function Base.get(pdb::ParameterDict, key, default)
+end
+
+
+function summarize(pdb::ParameterDict)
+ # compare pdb.provided, pdb.consumed, pdb.defaulted and display which
+ # parameters were unused, which ones were consumed and which ones
+ # falled back to defaults
+ # the output should also be logable to a file
+end
+
+
+# check if d.provided has key.
+# if it has one, copy it over to d.consumed and return it, if not copy the default
+# over to d.consumed and return it
+Base.get(d::ParameterDict, key, default) = get(d.provided, key, default)
+Base.length(d::ParameterDict) = length(d.provided)
+# similar to get but throw if key is not found in d.provided
+Base.getindex(d::ParameterDict, elements...) = getindex(d.dict, elements...)
+# disable this?
+Base.setindex!(d::ParameterDict, value, elements...) = setindex!(d.dict, value, elements...)
+# iteration interface
+Base.iterate(d::ParameterDict) = iterate(d.dict)
+Base.iterate(d::ParameterDict, index) = iterate(d.dict, index)
+# can this work with arbitrary dicts?
+Base.merge(d::ParameterDict, others::AbstractDict...) = ParameterDict(merge(d.dict, others...))
+# disable if we disable setindex!
+Base.pop!(d::ParameterDict, key, default) = pop!(d.dict, key, default)
+
+
+
+
+
+
+#######################################################################
+# main #
+#######################################################################
+
+
+parameters_db = ParameterDict(parameters)
+
+# factories, provided by dg1d, except equation should be provided by a project
+grid = make_Grid(parameters_db)
+cache = make_Cache(grid)
+hrsc = make_HRSC(parameters_db)
+equation = make_Equation(parameters_db)
+
+# we need an aggregation type
+disc = Discretization(grid, cache, hrsc, equation)
+
+# above might be combined into
+# and a project must then only implement make_Equation
+disc = Discretization(parameters_db)
+
+# prepare initial data, boundary conditions
+# this should also be implemented by a project
+# and just like in Makie, this should be merely are elaborated expression of what
+# problem we want to solve, e.g. which type of initial data we want and what boundary conditions
+# no actual computations should occur here
+problem = Problem(disc, parameters_db)
+
+# evolution
+runner = make_runner(disc, problem, parameters_db)
+
+summarize(runner, parameters_db)
+run(runner)
+
+function make_Runner(disc, problem, parameters_db)
+ # prepare evolution
+ rhs = make_rhs(disc, problem, parameters_db)
+ timestep = make_timestep(disc, problem, parameters_db)
+ callbacks = make_callbacks(disc, problem, parameters_db)
+ # () -> dg1d.evolve(rhs,
+end
+
+
+function make_rhs(disc, problem, parameters_db)
+ (u, t) -> dg1d.rhs(u, t, disc, problem, disc.equation)
+end
+
+
+function make_timestep(disc, problem, parameters_db)
+ (u, t) -> dg1d.timestep(u, t, disc, problem, disc.equation)
+end
+
+
+function dg1d.rhs(u, t, disc, problem, equation::ScalarEquation)
+ rhs(u, t, disc, problem, disc.hrsc)
+end
+
+
+function make_callbacks(disc, problem, parameters_db)
+ # setup rhs_callbacks, plot_callback, save_callback,
+ # id_smoother_callback, aligned_savecallback, restep_callback
+end
+
+
+# what needs to be implemented by a project will be
+# (I think we can't use factories because we cannot dispatch on any parameters then, or?)
+# - dg1d.make_equation
+# - dg1d.make_probelm
+# - dg1d.make_rhs
+# - dg1d.make_timestep
+# - dg1d.make_callback
+
+# things to figure out
+# - how to register variables consistently and in such a way that different
+# components can depend on one another? do we need another accepts/returns mechanism?
diff --git a/src/grid.jl b/src/grid.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c621ef08fddd21b7bb117eb03d09206cc5b8af95
--- /dev/null
+++ b/src/grid.jl
@@ -0,0 +1,281 @@
+#######################################################################
+# Type definition #
+#######################################################################
+
+
+struct Grid{N_Dim}
+ tree::Tree{N_Dim}
+ tensorbasis_cache::ObjectCache{Int64}
+ # interpolation_cache::ObjectCache{Int64}
+ bases::Vector{TensorBasis{N_Dim}}
+ ioffset::Vector{Int64}
+ ioffset_boundary::Vector{Int64}
+ boxes::Vector{Box{N_Dim}}
+ # parameters::Dict{String,Any}
+
+ function Grid{N_Dim}(ns, ps, extends, periodic, blueprint_basis) where N_Dim
+
+ tree = Tree{N_Dim}(ns..., periodic=periodic)
+ boxes = place_boxes(tree, extends...)
+
+ tensorbasis_cache = ObjectCache{Int64,SpectralElement}(blueprint_basis)
+ bases = [ TensorBasis([ tensorbasis_cache[p] for p in ps]...) for _ in 1:n_cells(tree) ]
+ ioffset, ioffset_boundary = Int64[1], Int64[1]
+ for b in bases
+ push!(ioffset, ioffset[end]+n_points(b))
+ push!(ioffset_boundary, ioffset_boundary[end]+n_points_boundary(b))
+ end
+
+ grid = new{N_Dim}(tree, tensorbasis_cache, bases, ioffset, ioffset_boundary, boxes)
+
+ return grid
+ end
+end
+const Grid1d = Grid{Cell1d}
+const Grid2d = Grid{Cell2d}
+const Grid3d = Grid{Cell3d}
+
+
+#######################################################################
+# Constructor #
+#######################################################################
+
+
+function Grid(parameters=Dict())
+ ns = query(Vector{Int64}, parameters, "ns", default=[5,5,5])
+ ps = query(Vector{Int64}, parameters, "ps", default=[ 5 for _ = 1:length(ns) ])
+ @assert all(ns .> 0)
+ @assert all(ps .== ps[1]) # equal polynomial orders in each direction for now
+ N_Dim, n_ps = length(ns), length(ps)
+ @assert N_Dim == n_ps """
+ Unable to determine dimension of grid. Please provide the same number of parameters for
+ - cells per dimension 'ns'
+ - polynomial orders per dimensions 'ps'
+ """
+
+ extends = query(Vector{Float64}, parameters, "extends", default=repeat([-1.0,1.0], N_Dim))
+ @assert length(extends) == 2*N_Dim """
+ Cannot determine domain size from 'extends' = $(extends).
+ Require one pair of extrema per dimension.
+ """
+ periodic = query(Vector{Bool}, parameters, "periodic", default=[ true for _ in 1:N_Dim ])
+ quadrature = query(String, parameters, "quadrature", default="LGL")
+
+ coordinates = query(String, parameters, "domaintype", default="cartesian")
+
+ blueprint(n) = SpectralElement(n, Symbol(quadrature))
+
+ return Grid{N_Dim}(ns, ps, (extends[2*(i-1)+1:2*(i-1)+2] for i = 1:N_Dim), tuple(periodic...), blueprint)
+end
+
+
+#######################################################################
+# Utilities #
+#######################################################################
+
+
+n_points(g::Grid) = sum(n_points.(g.tensorbasis))
+
+
+function Base.show(io::IO, grid::Grid)
+ ncs = ncells(grid)
+ npts = npoints(grid)
+ nbdrypts = npoints_boundary(grid)
+ print("""
+ $(typeof(grid))
+ # cells $(ncells(grid))
+ # points $(npoints(grid))
+ # boundary points $(npoints_boundary(grid))
+ """)
+end
+
+
+
+function sample_coordinates(grid::Grid{1})
+ xs = Float64[]
+ Js = Float64[]
+ detJ = Float64[]
+ npts = n_points(grid)
+ sizehint!(xs, npts)
+ sizehint!(Js, npts)
+ sizehint!(detJ, npts)
+ for (box, basis) in zip(grid.boxes, grid.tensorbasis)
+ (xmin, xmax), = extrema(basis)
+ map_x(z) = (xmax-xmin)/2 * z + (xmax+xmin)/2
+ box_xs = [ map_x(zi) for z in basis[Cart1d.X].z ]
+ box_Js = (xmax-xmin)/2 * ones(n_points(basis))
+ box_detJ = (xmax-xmin)/2
+ append!(xs, box_xs)
+ append!(Js, box_Js)
+ append!(detJ, box_detJ)
+ end
+ return [ xs ], [ Js ], detJ
+end
+
+
+function sample_coordinates(grid::Grid{2})
+ xs, ys = Float64[], Float64[]
+ Js_Xx, Js_Xy, Js_Yx, Js_Yy = [ Float64[] for _ = 1:4 ]
+ detJ = Float64
+ npts = n_points(grid)
+ sizehint!(xs, npts)
+ sizehint!(ys, npts)
+ sizehint!(Js_Xx, npts)
+ sizehint!(Js_Xy, npts)
+ sizehint!(Js_Yx, npts)
+ sizehint!(Js_Yy, npts)
+ sizehint!(detJ, npts)
+ for (box, basis) in zip(grid.boxes, grid.tensorbasis)
+ npts = n_points(basis)
+ (xmin, xmax), (ymin, ymax) = extrema(basis)
+ map_x(z) = (xmax-xmin)/2 * z + (xmax+xmin)/2
+ map_y(z) = (ymax-ymin)/2 * z + (ymax+ymin)/2
+ box_xs = [ map_x(zi) for z in basis[Cart2d.X].z ]
+ box_ys = [ map_y(zi) for z in basis[Cart2d.Y].z ]
+ box_xs, box_ys = mesh_coordinates(box_xs, box_ys)
+ box_Js_Xx = (xmax-xmin)/2 * ones(npts)
+ box_Js_Xy = zeros(npts)
+ box_Js_Yx = zeros(npts)
+ box_Js_Yy = (ymax-ymin)/2 * ones(npts)
+ box_detJ = @. box_Js_Xx * box_Js_Yy - box_Js_Xy * box_Js_Yx
+ push!(xs, box_xs)
+ push!(ys, box_ys)
+ push!(Js_Xx, box_Js_Xx)
+ push!(Js_Xy, box_Js_Xy)
+ push!(Js_Yx, box_Js_Yx)
+ push!(Js_Yy, box_Js_Yy)
+ push!(detJ, box_detJ)
+ end
+ return [ xs, ys ], [ Js_Xx, Js_Xy, Js_Yx, Js_Yy ], detJ
+end
+
+
+function sample_coordinates(grid::Grid{3})
+ xs, ys, zs = Float64[], Float64[], Float64[]
+ Js_Xx, Js_Xy, Js_Xz, Js_Yx, Js_Yy, Js_Yz, Js_Zx, Js_Zy, Js_Zz = [ Float64[] for _ = 1:9 ]
+ detJ = Float64[]
+ npts = n_points(grid)
+ sizehint!(xs, npts)
+ sizehint!(ys, npts)
+ sizehint!(zs, npts)
+ sizehint!(Js_Xx, npts)
+ sizehint!(Js_Xy, npts)
+ sizehint!(Js_Xz, npts)
+ sizehint!(Js_Yx, npts)
+ sizehint!(Js_Yy, npts)
+ sizehint!(Js_Yz, npts)
+ sizehint!(Js_Zx, npts)
+ sizehint!(Js_Zy, npts)
+ sizehint!(Js_Zz, npts)
+ sizehint!(detJ, npts)
+ for (box, basis) in zip(grid.boxes, grid.tensorbasis)
+ npts = n_points(basis)
+ (xmin, xmax), (ymin, ymax), (zmin, zmax) = extrema(basis)
+ map_x(z) = (xmax-xmin)/2 * z + (xmax+xmin)/2
+ map_y(z) = (ymax-ymin)/2 * z + (ymax+ymin)/2
+ map_z(z) = (zmax-zmin)/2 * z + (zmax+zmin)/2
+ box_xs = [ map_x(zi) for z in basis[Cart3d.X].z ]
+ box_ys = [ map_y(zi) for z in basis[Cart3d.Y].z ]
+ box_zs = [ map_z(zi) for z in basis[Cart3d.Z].z ]
+ box_xs, box_ys, box_zs = mesh(box_xs, box_ys, box_zs)
+ box_Js_Xx = (xmax-xmin)/2 * ones(npts)
+ box_Js_Xy = zeros(npts)
+ box_Js_Xz = zeros(npts)
+ box_Js_Yx = zeros(npts)
+ box_Js_Yy = (ymax-ymin)/2 * ones(npts)
+ box_Js_Yz = zeros(npts)
+ box_Js_Zx = zeros(npts)
+ box_Js_Zy = zeros(npts)
+ box_Js_Zz = (zmax-zmin)/2 * ones(npts)
+ box_detJ = ( @. box_Js_Xx * box_Js_Yy * box_Zz + box_Js_Xy * box_Js_Yz * box_Js_Zx +
+ box_Js_Xz * box_Js_Yx * box_Zy - box_Js_Xz * box_Js_Yy * box_Js_Zx -
+ box_Js_Xy * box_Js_Yx * box_Zz - box_Js_Xx * box_Js_Yz * box_Js_Zy )
+ append!(xs, box_xs)
+ append!(ys, box_ys)
+ append!(zs, box_zs)
+ push!(Js_Xx, box_Js_Xx)
+ push!(Js_Xy, box_Js_Xy)
+ push!(Js_Xz, box_Js_Xz)
+ push!(Js_Yx, box_Js_Yx)
+ push!(Js_Yy, box_Js_Yy)
+ push!(Js_Yz, box_Js_Yz)
+ push!(Js_Zx, box_Js_Zx)
+ push!(Js_Zy, box_Js_Zy)
+ push!(Js_Zz, box_Js_Zz)
+ push!(detJ, box_detJ)
+ end
+ return [ xs, ys, zs ], [ Js_Xx, Js_Xy, Js_Xz, Js_Yx, Js_Yy, Js_Yz, Js_Zx, Js_Zy, Js_Zz ], detJ
+end
+
+
+function mesh_coordinates(x, y)
+ nx, ny = length(x), length(y)
+ xx, yy = [ zeros(nx*ny) for _ in 1:2 ]
+ for i in 1:nx, j in 1:ny
+ ij = i + (j-1) * nx
+ xx[ij] = x[i]
+ yy[ij] = y[j]
+ end
+ return xx, yy
+end
+
+
+function mesh_coordinates(x,y,z)
+ nx, ny, nz = length(x), length(y), length(z)
+ xx, yy, zz = [ zeros(nx*ny*nz) for _ in 1:3 ]
+ for i in 1:nx, j in 1:ny, k in 1:nz
+ ijk = i + (j-1) * nx + (k-1) * nx * ny
+ xx[ijk] = x[i]
+ yy[ijk] = y[j]
+ zz[ijk] = z[k]
+ end
+ return xx, yy, zz
+end
+
+
+function register_cartesian_coordinates!(grid::Grid)
+ xys, components_Js, detJ = sample_coordinates(grid)
+ add_variables!(grid.cache, static_variables = [ xs..., Js..., detJ ])
+end
+
+
+function regrid!(grid::Grid)
+ # duplicate tree, bases and boxes and update them
+ # regrid cached variables inplace using the above copies
+ # e.g. perform various interpolation coarsening steps in different types of variables
+end
+
+
+function eachcell(grid::Grid, vars::RefBlockVector...)
+ ((block(v, i) for i in 1:n_cells(grid)) for v in vars)
+end
+
+
+# function range(grid::Grid, i)
+# range(grid.ioffset[i], length=n_points(grid.bases[i]))
+# end
+
+
+# struct Boundary{T<:Grid{N_Dim}}
+# grid::T
+# dir::Cart{N_Dim}.Direction
+# end
+#
+#
+# function eachcell(g::Boundary{Grid}, vars...)
+# ((view(v, range(grid, i, dir)) for v in vars) for i in 1:n_cells(grid))
+# end
+#
+#
+# function range(cell::Grid, i, dir)
+# range(bdry.grid.ioffset[i]+offset_boundary(grid.bases[i], dir),
+# length=n_points_boundary(grid.bases[i], dir))
+# end
+
+
+# for (ai,bi,ci) in eachcell(grid,a,b,c)
+# ci .= ai .+ bi
+# end
+# for (ai,bi,ci) in eachcell(Boundary(grid,Cart2d.Xmin),a,b,c)
+# ci .= ai .+ bi
+# end
diff --git a/src/grid.jl.old b/src/grid.jl.old
new file mode 100644
index 0000000000000000000000000000000000000000..5a79761715d3f36ed259138a6e76b4e4644e635a
--- /dev/null
+++ b/src/grid.jl.old
@@ -0,0 +1,293 @@
+@enum CellType::UInt8 VolumeCell BoundaryCell NoCell
+
+@enum Directions1d::UInt8 Left=1 Right
+@enum Directions2d::UInt8 Left=1 Right Top Bottom
+@enum Directions3d::UInt8 Left=1 Right Top Bottom Front Back
+
+Base.to_index(d::Directions1d) = Int(d)
+Base.to_index(d::Directions2d) = Int(d)
+Base.to_index(d::Directions3d) = Int(d)
+
+struct Cell{N_Dim,T_Neighbours}
+ type::CellType
+ index::Int64
+ level::Int64
+ neighbours::T_Neighbours
+ # hlevel::Int64
+ # plevel::Int64
+ # tlevel::Int64
+end
+
+const Cell1d = Cell{1,Neighbours1d}
+const Cell2d = Cell{2,Neighbours2d}
+const Cell3d = Cell{3,Neighbours3d}
+
+# struct CellId
+# index::Int64
+# type::CellType
+# end
+
+
+struct Tile{N_Dim,CellId}
+ cellids::NTuple{N_Dim,CellId}
+end
+const Tile1d{CellId} = Tile{1,CellId}
+const Tile2d{CellId} = Tile{2,CellId}
+
+
+is_boundary_cell(cid::CellId) = cid.type == BoundaryCell
+is_volume_cell(cid::CellId) = cid.type == VolumeCell
+is_no_cell(cid::CellId) = cid.type == NoCell
+is_boundary_cell(tcid::Tile{CellId}) = any(Tuple(cid.type == BoundaryCell for cid in tcid.cellids))
+is_volume_cell(tcid::Tile{CellId}) = any(Tuple(cid.type == VolumeCell for cid in tcid.cellids))
+is_no_cell(tcid::Tile{CellId}) = any(Tuple(cid.type == NoCell for cid in tcid.cellids))
+
+
+Base.length(cid::CellId) = 1
+Base.length(tcid::Tile1d) = 2
+Base.length(tcid::Tile2d) = 4
+
+
+const Neighbours1d = MVector{2,CellId}
+const Neighbours2d = MVector{4,Union{CellId,Tile1d}}
+const Neighbours3d = MVector{6,Union{CellId,Tile2d}}
+
+
+# # TODO Why not FieldVector? That would also allow indexing it.
+# mutable struct Neighbours1d
+# left::CellId
+# right::CellId
+# end
+#
+#
+# mutable struct Neighbours2d
+# left::Union{CellId,Tile{CellId}}
+# right::Union{CellId,Tile{CellId}}
+# top::Union{CellId,Tile{CellId}}
+# bottom::Union{CellId,Tile{CellId}}
+# end
+#
+#
+# mutable struct Neighbours3d
+# left::Union{CellId,Tile{CellId}}
+# right::Union{CellId,Tile{CellId}}
+# top::Union{CellId,Tile{CellId}}
+# bottom::Union{CellId,Tile{CellId}}
+# front::Union{CellId,Tile{CellId}}
+# back::Union{CellId,Tile{CellId}}
+# end
+
+
+# h-refinement
+has_tiles(nghbr::Neighbours1d) = false
+# has_tilled(nghbr::Neighbours) = any( #= member isa Tile for memeber in nghbr =#)
+find_tiles(nghbr::Neighbours1d) = nothing
+# findtilled(nghbr::Neighbours) = [ (side, id) for (side, id) in nghbr ]
+
+
+mutable struct Basis{N_Dim}
+ # TODO MVector instead of NTuple avoids making Basis mutable
+ elements::NTuple{N_Dim,SpectralElement}
+end
+Base.ndims(::Type{Basis{N_Dim}}) where N_Dim = N_Dim
+Base.ndims(b::Basis) = ndims(typeof(b))
+Base.size(basis::Basis) = Tuple( el.Npts for el in basis.elements )
+npoints(basis::Basis) = prod(size(basis))
+npoints_boundary(basis::Basis{1}) = 2
+function npoints_boundary(basis::Basis{2})
+ sz = size(basis)
+ return sum(2 .* sz) - 4 # remove double counted corners
+end
+function npoints_boundary(basis::Basis{3})
+ sz = size(basis)
+ return sum(2 .* (sz[1]*sz[2], sz[1]*sz[3], sz[2]*sz[3])) - 16 # remove double counted corners
+end
+
+
+mutable struct Cell{N_Dim,T_Neighbours}
+ # TODO replace (this,offset,isboundary) with CellId?
+ this::Int64
+ offset::Int64
+ neighbours::T_Neighbours
+ isboundary::Bool
+ basis::Basis{N_Dim}
+
+ # TODO incomplete constructor really needed?
+ Cell{N_Dim,T_Neighbours}() where {N_Dim,T_Neighbours} = new{N_Dim,T_Neighbours}()
+end
+Base.ndims(::Type{Cell{N_Dim,T_Neighbours}}) where {N_Dim,T_Neighbours} = N_Dim
+npoints(cell::Cell) = npoints(cell.basis)
+npoints_boundary(cell::Cell) = npoints_boundary(cell.basis)
+nneighbours(::Cell{N_Dim}) where N_Dim = 2*N_Dim
+
+
+const Cell1d = Cell{1,Neighbours1d}
+const Cell2d = Cell{2,Neighbours2d}
+const Cell3d = Cell{3,Neighbours3d}
+
+
+indices(cell::Cell) = range(cell.offset; length=npoints(cell))
+
+
+struct BoundaryIndices{N_Dim,T_Cell}
+ cell::T_Cell
+end
+const BoundaryIndices1d = BoundaryIndices{ndims(Cell1d),Cell1d}
+const BoundaryIndices2d = BoundaryIndices{ndims(Cell2d),Cell2d}
+const BoundaryIndices3d = BoundaryIndices{ndims(Cell3d),Cell3d}
+function boundary_indices(cell::Cell1d, side::Symbol)
+ @unpack offset = cell
+ if side === :left
+ return range(offset; length=1)
+ elseif side === :right
+ return range(offset+npoints(cell); length=1)
+ else
+ error()
+ end
+end
+function boundary_indices(cell::Cell2d, side::Symbol)
+ @unpack offset = cell
+ sz = size(cell)
+ if side === :left
+ return range(offset; length=sz[1])
+ elseif side === :right
+ return range(offset+npoints(cell)-sz[1]; length=sz[1])
+ elseif side === :top
+ return range(offset; length=sz[2], step=sz[1])
+ elseif side === :bottom
+ return range(offset+sz[1]; length=sz[2], step=sz[1])
+ else
+ error()
+ end
+end
+function boundary_indices(cell::Cell3d, side::Symbol)
+ @unpack offset = cell
+ sz = size(cell)
+ if side === :left
+ return range(offset; length=sz[1]*sz[2], step=1)
+ elseif side === :right
+ return range(offset+npoints(cell)-sz[1]*sz[2]; length=sz[1]*sz[2])
+ elseif side === :top
+ return range(offset; length=sz[1]*sz[2], step=sz[1])
+ elseif side === :bottom
+ return range(offset+sz[1]; length=sz[1]*sz[2], step=sz[1])
+ elseif side === :front
+ return range(offset+sz[1]*sz[2]; length=sz[2]*sz[3], step=sz[1]*sz[2])
+ elseif side === :back
+ return range(offset+npoints(cell)-sz[1]*sz[2]; length=sz[2]*sz[3], step=sz[1]*sz[2])
+ else
+ error()
+ end
+end
+
+
+struct Grid{T_Cell}
+ cells::Vector{T_Cell}
+end
+
+
+const Grid1d = Grid{Cell1d}
+const Grid2d = Grid{Cell2d}
+const Grid3d = Grid{Cell3d}
+
+
+function Grid(ncells)
+ grid = Grid1d([ Cell1d() for _ in 1:ncells ])
+ @unpack cells = grid
+ # connect cells
+ for (i, cell) in enumerate(cells)
+ i_left, i_right = periodize(i-1, ncells), periodize(i+1, ncells)
+ neighbours = Neighbours1d( CellId(i_left, VolumeCell), CellId(i_right, VolumeCell) )
+ setproperty!(cell, :neighbours, neighbours)
+ setproperty!(cell, :this, i)
+ end
+ return grid
+end
+
+
+function Grid(ncells_x, ncells_y)
+ ncells = ncells_x * ncells_y
+ grid = Grid2d([ Cell2d() for _ in 1:ncells ])
+ @unpack cells = grid
+ CI = CartesianIndices((ncells_x, ncells_y))
+ for (i, cell) in enumerate(cells)
+ ix, iy = Tuple(CI[i])
+ i_left, i_right = periodize(ix-1, ncells_x), periodize(ix+1, ncells_x)
+ i_top, i_bottom = periodize(iy-1, ncells_y), periodize(iy+1, ncells_y)
+ neighbours = Neighbours2d( CellId(i_left, VolumeCell), CellId(i_right, VolumeCell),
+ CellId(i_top, VolumeCell), CellId(i_bottom, VolumeCell) )
+ setproperty!(cell, :neighbours, neighbours)
+ setproperty!(cell, :this, i)
+ end
+ return grid
+end
+
+
+function Grid(ncells_x, ncells_y, ncells_z)
+ ncells = ncells_x * ncells_y * ncells_z
+ grid = Grid3d([ Cell3d() for _ in 1:ncells ])
+ @unpack cells = grid
+ CI = CartesianIndices((ncells_x, ncells_y, ncells_z))
+ for (i, cell) in enumerate(cells)
+ ix, iy, iz = Tuple(CI[i])
+ i_left, i_right = periodize(ix-1, ncells_x), periodize(ix+1, ncells_x)
+ i_top, i_bottom = periodize(iy-1, ncells_y), periodize(iy+1, ncells_y)
+ i_front, i_back = periodize(iz-1, ncells_z), periodize(iz+1, ncells_z)
+ neighbours = Neighbours3d( CellId(i_left, VolumeCell), CellId(i_right, VolumeCell),
+ CellId(i_top, VolumeCell), CellId(i_bottom, VolumeCell),
+ CellId(i_front, VolumeCell), CellId(i_back, VolumeCell) )
+ setproperty!(cell, :neighbours, neighbours)
+ setproperty!(cell, :this, i)
+ end
+ return grid
+end
+
+
+cells(grid::Grid) = grid.cells
+Base.ndims(grid::Grid) = ndims(eltype(cells(grid)))
+ncells(grid::Grid) = length(cells(grid))
+npoints(grid::Grid) = sum(npoints(c) for c in cells(grid))
+npoints_boundary(grid::Grid1d) = 2
+npoints_boundary(grid::Grid) = sum(npoints_boundary(c) for c in cells(grid))
+nneighbours(grid::Grid) = sum(nneighbours(c) for c in cells(grid))
+
+
+const _CACHE_SpectralElements = Dict{Int64,SpectralElement}()
+
+
+function make_SpectralElement(; args...)
+ N = args[:N]
+ if haskey(_CACHE_SpectralElements, N)
+ return _CACHE_SpectralElements[N]
+ end
+ el = SpectralElement(values(args)...)
+ _CACHE_SpectralElements[N] = el
+ return el
+end
+
+
+function Base.fill!(grid::Grid, basis::Basis)
+ @assert ndims(basis) == ndims(grid)
+ npts = npoints(basis)
+ offset = 1
+ for cell in cells(grid)
+ setproperty!(cell, :basis, basis)
+ setproperty!(cell, :offset, offset)
+ offset += npts
+ end
+end
+
+
+function Base.show(io::IO, grid::Grid)
+ ncs = ncells(grid)
+ npts = npoints(grid)
+ nbdrypts = npoints_boundary(grid)
+ print("""
+ $(typeof(grid))
+ # cells $(ncells(grid))
+ # points $(npoints(grid))
+ # boundary points $(npoints_boundary(grid))
+ """)
+end
+
+# TODO Show methods for cell and basis.
diff --git a/src/ideas/cache.jl b/src/ideas/cache.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d8ba1dbb307dafe1cec94a9e8c36a7f7717b7137
--- /dev/null
+++ b/src/ideas/cache.jl
@@ -0,0 +1,144 @@
+export VariableCache, register_variables!, add_variables!, wrap_variables!
+
+
+#######################################################################
+# Types #
+#######################################################################
+
+
+"""
+ VariableVariableCache
+
+Caching object to bundle variables of different kinds.
+"""
+struct VariableCache{T_Grid<:AbstractGrid}
+ grid::T_Grid
+ varlist::Dict{Symbol,Symbol}
+ dynamic::Variables
+ static::Variables
+ cell::Variables
+ rhs::Variables
+ # interface::Variables # should we added this to separate inferface and boundary values
+ # would allow to easily distinguish boundary values that are computed from boundary conditions
+ # and/or neighboring grids (e.g. sphered cubes etc).
+ boundary::Variables
+ domain::Variables
+end
+
+
+function VariableCache(grid::Grid)
+ verbose = false
+ varlist = Dict{Symbol, Symbol}()
+ npts, npts_bdry = n_points(grid), n_points_boundary(grid)
+ ncells = n_cells(grid)
+ dynamic = Variables(npts, verbose=verbose)
+ static = Variables(npts, verbose=verbose)
+ cell = Variables(ncells, verbose=verbose)
+ rhs = Variables(npts, verbose=verbose)
+ boundary = Variables(npts_bdry, verbose=verbose)
+ domain = Variables(1, verbose=verbose)
+ return VariableCache(grid, varlist, dynamic, static, cell, rhs, boundary, domain)
+end
+
+
+#######################################################################
+# Methods #
+#######################################################################
+
+
+@inline UnPack.unpack(x::VariableCache, ::Val{f}) where {f} = get_variable(x, Val(f))
+
+
+function get_variable(varcache::VariableCache, ::Val{f}) where {f}
+ field = v.varlist[f]
+ return getproperty(varcache, field)
+end
+
+
+
+
+# TODO Can we avoid this and instead overload Unpack.@unpack or so?
+function Base.propertynames(::VariableCache)
+ fnames = fieldnames(VariableCache)
+ filter(f -> !in(f, [ :varlist ]), fnames)
+end
+
+
+"""
+ wrap_variables!(C::VariableCache, u=nothing, du=nothing)
+
+Wrap fields `dynamic, rhs` around data `u, du` that is not owned by `C`.
+"""
+function wrap_variables!(C::VariableCache; u=nothing, du=nothing)
+ !isnothing(u) && set!(C.dynamic, u)
+ !isnothing(du) && set!(C.rhs, du)
+ return
+end
+
+
+
+function register_variables!(C::VariableCache;
+ dynamic = [], rhs = [], static = [], cell = [], bdry = [], domain = [])
+
+ register!(C.dynamic, Symbol.(dynamic)...)
+ register!(C.rhs, Symbol.(rhs)...)
+ register!(C.static, Symbol.(static)...)
+ register!(C.cell, Symbol.(cell)...)
+ register!(C.boundary, Symbol.(boundary)...)
+ register!(C.domain, Symbol.(domain)...)
+
+ empty!(cache.varlist)
+ push_names_field!(cache.varlist, dg1d.names(cache.dynamic), :dynamic)
+ push_names_field!(cache.varlist, dg1d.names(cache.rhs), :rhs)
+ push_names_field!(cache.varlist, dg1d.names(cache.static), :static)
+ push_names_field!(cache.varlist, dg1d.names(cache.cell), :cell)
+ push_names_field!(cache.varlist, dg1d.names(cache.boundary), :boundary)
+ push_names_field!(cache.varlist, dg1d.names(cache.domain), :domain)
+
+ return
+end
+
+
+push_names_field!(dict, names, field) = for name in names; push!(dict, name => field); end
+
+
+function Base.show(io::IO, C::VariableCache)
+ variables_list = collect(keys(C.fields_to_variables))
+ N_variables = length(variables_list)
+ if N_variables == 0
+ @warn "No variables registered!"
+ return
+ end
+ names_list = collect(values(C.fields_to_variables))
+ permute_idx = sortperm(variables_list)
+ variables_list, names_list = variables_list[permute_idx], names_list[permute_idx]
+ max_rows = maximum(length.(names_list))
+ data = Matrix{String}(undef, max_rows, N_variables)
+ fill!(data, "")
+ for (idx, names) in enumerate(names_list)
+ N_names = length(names)
+ data[1:N_names,idx] .= string.(names)
+ end
+ # sort!(data, dims=1)
+ pretty_table(data, header=variables_list, crop=:none)
+end
+
+
+function resize!(varcache::VariableCache)
+ npts = n_points(varcache)
+ npts_bdry = n_points_boundary(varcache)
+ ncells = n_cells(varcache)
+ resize!(varcache.dynamic)
+ resize!(varcache.rhs, npts)
+ resize!(varcache.static, npts)
+ resize!(varcache.cell, ncells)
+ resize!(varcache.boundary, npts_bdry)
+ # domain unchanged, because it only contains single values
+end
+
+
+function regrid_dynamic!(varcache::VariableCache)
+ for name in varcache.dynamic
+
+ end
+end
diff --git a/src/ideas/grhd_rhs.jl b/src/ideas/grhd_rhs.jl
new file mode 100644
index 0000000000000000000000000000000000000000..3ddf94219967e676fd78753e0f05e6531bb26c57
--- /dev/null
+++ b/src/ideas/grhd_rhs.jl
@@ -0,0 +1,136 @@
+struct Grid{N_Dim,T_SpectralElement}
+ mesh::TreeMesh{N_Dim}
+ bases_cache::ObjectCache{Int64,T_SpectralElement}
+ bases::Vector{TensorBasis{N_Dim}}
+ cache::Cache
+end
+
+
+@unpack variable"u, v, w, fu, fv, fw" = cache
+@unpack dynamic"u, v, w" = cache
+@unpack static"u, v, w" = cache
+@unpack boundary"u, v, w" = cache
+@unpack cell"u, v, w" = cache
+# register_variable.(cache, dynamic"u, v, w")
+->
+getproperty(getproperty(cache, :dynamic), :u)
+(alternative) ->
+getproperty(cache, Dynamic(:u))
+(intermediate) ->
+@unpack Var(:u), Var(:v), Var(:w), ... = cache
+
+
+function Grid(tm::TreeMesh{N_Dim},
+ grid_extends::Vararg{Tuple{Float64,Float64},N_Dim}) where N_Dim
+
+ ncells = n_cells(tm)
+ bases = [ TensorBasis(plevels(c)) for c in cells(tm) ]
+end
+
+
+
+const Grid1d = Grid{Cell1d}
+const Grid2d = Grid{Cell2d}
+const Grid3d = Grid{Cell3d}
+
+
+function Grid(ncells)
+ grid = Grid1d([ Cell1d() for _ in 1:ncells ])
+ @unpack cells = grid
+ # connect cells
+ for (i, cell) in enumerate(cells)
+ i_left, i_right = periodize(i-1, ncells), periodize(i+1, ncells)
+ neighbours = Neighbours1d( CellId(i_left, VolumeCell), CellId(i_right, VolumeCell) )
+ setproperty!(cell, :neighbours, neighbours)
+ setproperty!(cell, :this, i)
+ end
+ return grid
+end
+
+
+function Grid(ncells_x, ncells_y)
+ ncells = ncells_x * ncells_y
+ grid = Grid2d([ Cell2d() for _ in 1:ncells ])
+ @unpack cells = grid
+ CI = CartesianIndices((ncells_x, ncells_y))
+ for (i, cell) in enumerate(cells)
+ ix, iy = Tuple(CI[i])
+ i_left, i_right = periodize(ix-1, ncells_x), periodize(ix+1, ncells_x)
+ i_top, i_bottom = periodize(iy-1, ncells_y), periodize(iy+1, ncells_y)
+ neighbours = Neighbours2d( CellId(i_left, VolumeCell), CellId(i_right, VolumeCell),
+ CellId(i_top, VolumeCell), CellId(i_bottom, VolumeCell) )
+ setproperty!(cell, :neighbours, neighbours)
+ setproperty!(cell, :this, i)
+ end
+ return grid
+end
+
+
+function Grid(ncells_x, ncells_y, ncells_z)
+ ncells = ncells_x * ncells_y * ncells_z
+ grid = Grid3d([ Cell3d() for _ in 1:ncells ])
+ @unpack cells = grid
+ CI = CartesianIndices((ncells_x, ncells_y, ncells_z))
+ for (i, cell) in enumerate(cells)
+ ix, iy, iz = Tuple(CI[i])
+ i_left, i_right = periodize(ix-1, ncells_x), periodize(ix+1, ncells_x)
+ i_top, i_bottom = periodize(iy-1, ncells_y), periodize(iy+1, ncells_y)
+ i_front, i_back = periodize(iz-1, ncells_z), periodize(iz+1, ncells_z)
+ neighbours = Neighbours3d( CellId(i_left, VolumeCell), CellId(i_right, VolumeCell),
+ CellId(i_top, VolumeCell), CellId(i_bottom, VolumeCell),
+ CellId(i_front, VolumeCell), CellId(i_back, VolumeCell) )
+ setproperty!(cell, :neighbours, neighbours)
+ setproperty!(cell, :this, i)
+ end
+ return grid
+end
+
+
+cells(grid::Grid) = grid.cells
+Base.ndims(grid::Grid) = ndims(eltype(cells(grid)))
+ncells(grid::Grid) = length(cells(grid))
+npoints(grid::Grid) = sum(npoints(c) for c in cells(grid))
+npoints_boundary(grid::Grid1d) = 2
+npoints_boundary(grid::Grid) = sum(npoints_boundary(c) for c in cells(grid))
+nneighbours(grid::Grid) = sum(nneighbours(c) for c in cells(grid))
+
+
+const _CACHE_SpectralElements = Dict{Int64,SpectralElement}()
+
+
+function make_SpectralElement(; args...)
+ N = args[:N]
+ if haskey(_CACHE_SpectralElements, N)
+ return _CACHE_SpectralElements[N]
+ end
+ el = SpectralElement(values(args)...)
+ _CACHE_SpectralElements[N] = el
+ return el
+end
+
+
+function Base.fill!(grid::Grid, basis::Basis)
+ @assert ndims(basis) == ndims(grid)
+ npts = npoints(basis)
+ offset = 1
+ for cell in cells(grid)
+ setproperty!(cell, :basis, basis)
+ setproperty!(cell, :offset, offset)
+ offset += npts
+ end
+end
+
+
+function Base.show(io::IO, grid::Grid)
+ ncs = ncells(grid)
+ npts = npoints(grid)
+ nbdrypts = npoints_boundary(grid)
+ print("""
+ $(typeof(grid))
+ # cells $(ncells(grid))
+ # points $(npoints(grid))
+ # boundary points $(npoints_boundary(grid))
+ """)
+end
+
+# TODO Show methods for cell and basis.
diff --git a/src/ideas/mesh2d.jl b/src/ideas/mesh2d.jl
new file mode 100644
index 0000000000000000000000000000000000000000..9773e0fd8e45bf708bec998b61b6ed0915b00b65
--- /dev/null
+++ b/src/ideas/mesh2d.jl
@@ -0,0 +1,372 @@
+export Mesh, make_Mesh
+
+
+#######################################################################
+# Types #
+#######################################################################
+
+
+abstract type AbstractMesh{Ndim} end
+
+@doc """
+ Mesh{N}
+
+Object that bundles grid and operators needed for modeling a Discontinuous Galerkin
+scheme in N dimensions.
+"""
+struct Mesh{Ndim} <: AbstractMesh{Ndim}
+ elements::NTuple{Ndim,SpectralElement} # collection of spectral element operators
+ Ks::NTuple{Ndim,Int64} # number of cells
+ bounds::NTuple{Ndim,NTuple{Ndim,Float64}}
+ CFL::Float64 # Courant Friedrichs Lewy factor
+ Ls::NTuple{Ndim,Float64} # domain length
+ dls::NTuple{Ndim,Float64} # equidistant subcell length
+ xs::NTuple{Ndim,Matrix{Float64}} # grid points
+ invjacs::NTuple{Ndim,Matrix{Float64}} # inverse jacobian of coordinate transformation
+ # from computational [-1,1] onto physical grid [xmin,xmax]
+
+ function Mesh{Ndim}(; Ns, Ks, bounds, CFL, quadrature_method)
+ elements = [ SpectralElement(N, Symbol(quadrature_method)) for N in Ns ]
+ zs = [ el.z for el in elements ]
+ bounds = [ extrema(el.z) for el in elements ]
+ Ls = [ diff(b) for b in bounds ]
+ dls = [ L/K for (L,K) in zip(Ls, Ks) ]
+ as = [ [ b[1] + dl * idx for idx = 0 : K-1 ] for (b,dl,K) in zip(bounds, Ls, Ks) ]
+ bs = [ a.+dl for (a,dl) in zip(as,dls) ]
+ Npts = [ length(z) for z in zs ]
+ xs = [ Matrix{Float64}(undef, Np, K) for (Np,K) in zip(Npts, Ks) ]
+ invjacs = [ similar(x) for x in xs ]
+ for (Np, K, a, b, z, x, invjac) in zip (Npts, Ks, as, bs, zs, xs, invjacs)
+ for i=1:Np, j=1:K
+ aj = a[j]
+ bj = b[j]
+ zi = z[i]
+ x[i,j] = 0.5 * ((bj - aj) * zi + (aj + bj))
+ invjac[i,j] = 2.0 / (bj - aj)
+ end
+ # sanitize cell end points to have same x coordinate at common interfaces
+ for j=1:K-1
+ xl = x[end,j]
+ xr = x[1,j+1]
+ xmean = 0.5 * (xl + xr)
+ x[end,j] = xmean
+ x[1,j+1] = xmean
+ end
+ end
+ return new{Ndim}(elements, Tuple(Ks), bounds, CFL, Ls, dls, xs, invjacs)
+ end
+end
+
+
+#######################################################################
+# Factories #
+#######################################################################
+
+
+make_Mesh(prms) = Mesh(; parse_parameters(Mesh, prms)...)
+make_Mesh(; args...) = make_Mesh(FlatDict(args))
+
+
+######################################################################
+# Parameters #
+######################################################################
+
+
+const Mesh1d = Mesh{1}
+const Mesh2d = Mesh{2}
+const Mesh3d = Mesh{3}
+
+
+function parse_parameters(::Type{Mesh{Ndim}}, prms)
+ N = query(Vector{Int64}, prms, "N", default=repeat([5],Ndim))
+ K = query(Vector{Int64}, prms, "K", default=repeat([10],Ndim))
+ CFL = query(Float64, prms, "CFL")
+ quadrature_method = query(String, prms, "quadrature_method", default="LGL",
+ options=["LGL", "GLGL", "LGL_lumped"])
+ range = query(Vector{Vector{Float64}}, prms, "range", default=repeat([[-1.0,1.0]], Ndim))
+ @assert all(N .> 0) "'N > 0' required"
+ @assert all(K .> 0) "'K > 0' required"
+ @assert CFL > 0 "'CFL > 0' required"
+ @assert all(length.(range) .== 2) "'range' must be a tuple"
+ # @assert all(range[1] < xrange[2] "'xrange' limits must be ordered"
+ @assert quadrature_method in ["LGL", "GLGL", "LGL_lumped"] "'quadrature_method' invalid"
+ return N, K, CFL, quadrature_method, range
+end
+
+
+function parse_parameters(::Type{Mesh1d}, prms)
+ N, K, CFL, quadrature_method, range = parse_parameters(Mesh, prms)
+
+
+end
+
+
+function parse_parameters(::Type{Mesh1d}, prms)
+ Ns = query(Int64, prms, "N", default=[5])
+ Ks = query(Int64, prms, "K", default=[20])
+ CFL = query(Float64, prms, "CFL", default=0.5)
+ quadrature_method = query(String, prms, "quadrature_method", default="LGL",
+ options=["LGL", "GLGL", "LGL_lumped"])
+ xrange = query(Vector{Float64}, prms, "xrange", default = [-1.0,1.0])
+ # sanitize
+ @assert Ns > 0 "'N > 0' required"
+ @assert Ks > 0 "'K > 0' required"
+ @assert CFL > 0 "'CFL > 0' required"
+ @assert length(xrange) == 2 "'xrange' must be a tuple"
+ @assert xrange[1] < xrange[2] "'xrange' limits must be ordered"
+ @assert quadrature_method in ["LGL", "GLGL", "LGL_lumped"] "'quadrature_method' invalid"
+ bounds = (xrange,)
+ return (; Ns, Ks, CFL, quadrature_method, bounds)
+end
+
+
+function parse_parameters(::Type{Mesh2d}, prms)
+ Ns = query(NTuple{2,Int64}, prms, "N", default=(5,5))
+ Ks = query(NTuple{2,Int64}, prms, "K", default=(20,20))
+ CFL = query(Float64, prms, "CFL", default=0.5)
+ quadrature_method = query(String, prms, "quadrature_method", default="LGL",
+ options=["LGL", "GLGL", "LGL_lumped"])
+ xrange = query(Vector{Float64}, prms, "xrange", default = [-1.0,1.0])
+ yrange = query(Vector{Float64}, prms, "yrange", default = [-1.0,1.0])
+ # sanitize
+ @assert all(Ns .> 0) "'N > 0' required"
+ @assert all(Ks .> 0) "'K > 0' required"
+ @assert CFL > 0 "'CFL > 0' required"
+ @assert length(xrange) == 2 "'xrange' must be a tuple"
+ @assert length(yrange) == 2 "'yrange' must be a tuple"
+ @assert xrange[1] < xrange[2] "'xrange' limits must be ordered"
+ @assert yrange[1] < yrange[2] "'yrange' limits must be ordered"
+ @assert quadrature_method in ["LGL", "GLGL", "LGL_lumped"] "'quadrature_method' invalid"
+ bounds = (xrange, yrange)
+ return (; Ns, Ks, CFL, quadrature_method, bounds)
+end
+
+
+function parse_parameters(::Type{Mesh3d}, prms)
+ Ns = query(NTuple{3,Int64}, prms, "N", default=(5,5,5))
+ Ks = query(NTuple{3,Int64}, prms, "K", default=(20,20,20))
+ CFL = query(Float64, prms, "CFL", default=0.5)
+ quadrature_method = query(String, prms, "quadrature_method", default="LGL",
+ options=["LGL", "GLGL", "LGL_lumped"])
+ xrange = query(Vector{Float64}, prms, "xrange", default = [-1.0,1.0])
+ yrange = query(Vector{Float64}, prms, "yrange", default = [-1.0,1.0])
+ zrange = query(Vector{Float64}, prms, "zrange", default = [-1.0,1.0])
+ # sanitize
+ @assert all(Ns .> 0) "'N > 0' required"
+ @assert all(Ks .> 0) "'K > 0' required"
+ @assert CFL > 0 "'CFL > 0' required"
+ @assert length(xrange) == 2 "'xrange' must be a tuple"
+ @assert length(yrange) == 2 "'yrange' must be a tuple"
+ @assert length(zrange) == 2 "'zrange' must be a tuple"
+ @assert xrange[1] < xrange[2] "'xrange' limits must be ordered"
+ @assert yrange[1] < yrange[2] "'yrange' limits must be ordered"
+ @assert zrange[1] < zrange[2] "'zrange' limits must be ordered"
+ @assert quadrature_method in ["LGL", "GLGL", "LGL_lumped"] "'quadrature_method' invalid"
+ bounds = (xrange, yrange, zrange)
+ return (; Ns, Ks, CFL, quadrature_method, bounds)
+end
+
+
+#######################################################################
+# Methods #
+#######################################################################
+
+
+Base.broadcastable(mesh::AbstractMesh) = Ref(mesh)
+Base.ndims(mesh::AbstractMesh{N}) = N
+coordinates(mesh::AbstractMesh) = mesh.xs
+bounds(mesh::AbstractMesh) = mesh.bounds
+Base.length(mesh::AbstractMesh{N})::NTuple{N,Int64} = prod(length.(coordinates(mesh)))
+n_points(mesh::AbstractMesh{N}) = mesh.element.Npts .+ 1
+n_cells(mesh::AbstractMesh{N}) = mesh.Ks
+layout(mesh::AbstractMesh{N})::NTuple{N,Tuple{Int,Int}} = (n_points(mesh), n_cells(mesh))
+
+
+"""
+ cell_bdrys(mesh::AbstarctMesh)
+
+Return coordinates of cell boundaries from a `mesh`.
+"""
+function cell_bdrys(mesh::Mesh{N}) where N
+ bdrys = NTuple{N,Vector{Float64}}([ Vector{Float64}(undef, K+1) for K in Ks ])
+ for (idx,xi) in enumerate(coordinates(mesh))
+ bdrys[idx] .= vcat(x[1,:], [x[end,end]])
+ end
+ return bdrys
+end
+
+@deprecate extract_cell_bdrys(mesh) cell_bdrys(mesh)
+
+@doc """
+ grid_average(u, mesh::Mesh)
+
+Return the grid average of all cell averages
+"""
+function grid_average(u, mesh::Mesh1d)
+ @unpack element, dl, L = mesh
+ cellavg = cellwise_average(u, mesh) .* dl
+ avg = sum(cellavg) / L
+ return avg
+end
+
+
+@doc """
+ `cellwise_average(u, mesh::Mesh)`
+
+Return the cellwise average of `u` over `mesh.x`.
+"""
+function cellwise_average(u, mesh::Mesh)
+ @unpack element, invjac, K, dl = mesh
+ @unpack Npts = element
+ v = ones(Npts)
+ avg = zeros(K)
+ for i in 1:K
+ avg[i] = element.iproduct(u[:,i], v)
+ end
+ avg ./= invjac[1,:] * dl # constant Jacobian, constant cell width
+ return avg
+end
+
+
+@doc """
+ `cellwise_inner_product(u, v, mesh::Mesh)`
+
+Return vector of inner products between `u, v` on every cell in `mesh`.
+"""
+function cellwise_inner_product(u, v, mesh::Mesh)
+ @unpack element, K, invjac = mesh
+ @unpack Npts = element
+ L = layout(mesh)
+ mat_u = reshape(view(u,:), L)
+ mat_v = reshape(view(v,:), L)
+ inner_products = [ element.iproduct(mat_u[:,k], mat_v[:,k]) for k = 1:K ]
+ inner_products ./= invjac[1,:] # constant Jacobian
+ return inner_products
+end
+
+
+@doc """
+ `broken_inner_product(u, v, mesh::Mesh)`
+
+Return sum of all inner products between `u,v` on every cell in `mesh`.
+"""
+function broken_inner_product(u, v, mesh::Mesh)
+ inner_products = cellwise_inner_product(u, v, mesh)
+ total = sum(inner_products[:])
+ return total
+end
+
+
+# function LinearAlgebra.norm(u, mesh::Mesh, order)
+# @match order begin
+# 1 => norm_L1(u, mesh)
+# 2 => norm_L2(u, mesh)
+# _ => error("Norm for order = $order not implemented.")
+# end
+# end
+# norm_L2(u, mesh::Mesh) = sqrt(broken_inner_product(u, u, mesh))
+# norm_L1(u, mesh::Mesh) = broken_inner_product(abs.(u), ones(Float64, size(u)), mesh)
+
+
+@doc """
+Compute local maximum absolute value.
+"""
+function local_maxabs(v::AbstractMatrix, mesh::Mesh)
+ return vec(maximum(abs.(v), dims=1))
+end
+
+
+@doc """
+Compute local minimum absolute value.
+"""
+function local_maxabs(v::AbstractVector, mesh::Mesh)
+ return local_maxabs(reshape(v, (mesh.element.Npts, mesh.K)))
+end
+
+
+"""
+ find_nearest(xi, mesh::Mesh1d)
+ find_nearest(xi, mesh::Mesh2d)
+ find_nearest(xi, mesh::Mesh3d)
+
+Find index of node in `mesh` that lies closest to `xi`.
+Returns the index of the outermost point if `xi` falls outside `mesh`.
+"""
+function find_nearest(pt_x, mesh::Mesh1d)
+ px = pt_x
+ x, = coordinates(mesh)
+ min_ix = 1
+ min_distance = abs(x[min_ix] - px)
+ for ix = 2:length(x)
+ distance = abs(x[ix] - px)
+ if distance < min_distance
+ min_distance = distance
+ min_ix = ix
+ end
+ end
+ return (min_ix,)
+end
+function find_nearest(pt_xy, mesh::Mesh2d)
+ px, py = pt_xy
+ xs, ys = coordinates(mesh)
+ min_ix, min_iy = 1, 1
+ min_distance = sqrt( (xs[min_ix]-px)^2 + (ys[min_iy]-py)^2 )
+ for ix = 2:length(x), iy = 2:length(y)
+ distance = sqrt( (xs[ix]-px)^2 + (ys[iy]-py)^2 )
+ if distance < min_distance
+ min_distance = distance
+ min_ix, min_iy = ix, iy
+ end
+ end
+ return (min_ix, min_iy)
+end
+function find_nearest(pt_xyz, mesh::Mesh2d)
+ px, py, pz = pt_xyz
+ xs, ys, zs = coordinates(mesh)
+ min_ix, min_iy, min_iz = 1, 1, 1
+ min_distance = sqrt( (xs[min_ix]-px)^2 + (ys[min_iy]-py)^2 + (zs[min_iz]-pz)^2 )
+ for ix = 2:length(x), iy = 2:length(y), iz = 2:length(iz)
+ distance = sqrt( (xs[ix]-px)^2 + (ys[iy]-py)^2 + (zs[min_iz]-pz)^2 )
+ if distance < min_distance
+ min_distance = distance
+ min_ix, min_iy, min_iz = ix, iy, iz
+ end
+ end
+ return (min_ix, min_iy, min_iz)
+end
+
+
+@deprecate locate_point(xi, mesh::Mesh) find_nearest(xi, mesh::Mesh)
+@deprecate locate_point(xi, mesh::AbstractMesh) find_nearest(xi, mesh::AbstractMesh)
+
+
+"""
+ find_cell_index(xi, mesh::Mesh)
+
+Find the cell index into which coordinate `xi` falls into the `mesh`.
+Note that it returns the first index found in case `xi` falls onto an interface boundary.
+Returns nothing if `xi` lies outside the `mesh`.
+"""
+function find_cell_index(pt_x, mesh::Mesh1d)
+ px = pt_x
+ xbdrys, = cell_bdrys(mesh)
+ return findfirst(1:length(xbdrys)-1) do ix
+ xbdrys[ix] <= px <= xbdrys[ix+1]
+ end
+end
+function find_cell_index(pt_xy, mesh::Mesh2d)
+ px, py = pt_xy
+ xbdrys, ybdrys = cell_bdrys(mesh)
+ return findfirst((1:length(xbdrys)-1,1:length(ybdrys)-1)) do ix, iy
+ xbdrys[ix] <= px <= xbdrys[ix+1] && ybdrys[iy] <= py <= ybdrys[iy+1]
+ end
+end
+function find_cell_index(pt_xyz, mesh::Mesh3d)
+ px, py, pz = pt_xyz
+ xbdrys, ybdrys, zbdrys = cell_bdrys(mesh)
+ return findfirst((1:length(xbdrys)-1,1:length(ybdrys)-1,1:length(zbdrys)-1)) do ix, iy, iz
+ xbdrys[ix] <= px <= xbdrys[ix+1] &&
+ ybdrys[iy] <= py <= ybdrys[iy+1] &&
+ zbdrys[iz] <= pz <= zbdrys[iz+1]
+ end
+end
diff --git a/src/ideas/srhd_rhs.jl b/src/ideas/srhd_rhs.jl
new file mode 100644
index 0000000000000000000000000000000000000000..a928f52de88952a2b1aa9e48b378d0666e8ee0a7
--- /dev/null
+++ b/src/ideas/srhd_rhs.jl
@@ -0,0 +1,321 @@
+Base.@kwdef struct RHS{T_Equation <:SREulerEquation,
+ T_RSolver <:AbstractRiemannSolver,
+ T_HRSC <:Maybe{HRSC.AbstractHRSC},
+ T_TCI <:Maybe{TCI.AbstractTCI},
+ T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
+ T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
+
+ cache::Cache
+ mesh::Mesh
+ equation::T_Equation
+ rsolver::T_RSolver
+ hrsc::T_HRSC = nothing
+ tci::T_TCI = nothing
+ ldg_rsolver::T_LDG_RSolver = nothing
+ av_rsolver::T_AV_RSolver = nothing
+end
+
+
+#######################################################################
+# Register variables #
+#######################################################################
+
+
+# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
+# definition in this module without a warning or error. This does not happen for methods
+# defined in Base.
+# Because of that we prefix with _ below.
+# Cf. https://github.com/JuliaLang/julia/issues/15510
+function _register_variables!(rhs::RHS)
+ register_variables!(rhs.cache,
+ dynamic_variablenames = (:D, :S, :tau,), # conservatives
+ rhs_variablenames = (:rhs_D, :rhs_S, :rhs_tau), # rhs of conservatives
+ static_variablenames = (:flx_D, :flx_S, :flx_tau, # conservatives' fluxes
+ :rho, :v, :p, :eps, # primitives
+ :max_v, :x, # maximum charactersitic speed, grid
+ :E, :Em1, # entropy
+ :flx_E, :flx_Em1), # entropy flux
+ cell_variablenames = (:cellmax_v,),
+ global_variablenames = (:t, :tm1) # time steps
+ )
+ _register_variables!(rhs, rhs.hrsc)
+ _register_variables!(rhs, rhs.tci)
+ @unpack x = get_static_variables(rhs.cache)
+ @. x = rhs.mesh.x[:]
+end
+
+
+function _register_variables!(rhs, tci_or_hrsc::Nothing) end
+
+
+_register_variables!(rhs, hrsc::AbstractHRSC) =
+ error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
+
+_register_variables!(rhs, hrsc::HRSC.AbstractReconstruction) = nothing
+
+
+function _register_variables!(rhs, hrsc::HRSC.AbstractArtificialViscosity)
+ register_variables!(rhs.cache,
+ static_variablenames = (:ldg_D, :ldg_S, :ldg_tau, # local DG variable
+ :flx_ldg_D, :flx_ldg_S, :flx_ldg_tau), # local DG flux
+ cell_variablenames = (:mu,), # 'one viscosity to rule them all'
+ )
+end
+
+
+function _register_variables!(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
+ _register_variables!(rhs, hrsc.av)
+ register_variables!(rhs.cache,
+ static_variablenames = (:smoothed_mu,)
+ )
+end
+
+
+function _register_variables!(rhs, tci::TCI.AbstractTCI)
+ register_variables!(rhs.cache, cell_variablenames = (:flag,:D_flag,:S_flag,:tau_flag))
+end
+
+
+function _register_variables!(rhs, tci::TCI.EntropyProduction)
+ register_variables!(rhs.cache,
+ cell_variablenames = (:flag,),
+ static_variablenames = (:EP,), # entropy production
+ bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
+ # entropy production due to jumps at interfaces
+ )
+end
+
+
+#######################################################################
+# Timestep #
+#######################################################################
+
+
+function timestep(state_u, t, rhs::RHS)
+ wrap_dynamic_variables!(rhs.cache, state_u)
+ return timestep(rhs, rhs.hrsc)
+end
+
+
+function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
+ @unpack cache, mesh, equation = rhs
+ @unpack max_v = get_static_variables(cache)
+
+ broadcast_volume!(cons2prim, equation, cache)
+ broadcast_volume!(speed, equation, cache)
+ vmax = dg1d.absolute_maximum(max_v)
+ # vmax > 0.99 && error("Unreasonable maximum speed encountered: $vmax")
+ vmax_limit = 0.9999
+ if vmax > vmax_limit
+ @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
+ vmax = vmax_limit
+ end
+
+ @unpack dl, x, CFL, element = mesh
+ @unpack N = element
+
+ dt = CFL * dl / (N^2 * vmax)
+ return dt
+end
+
+
+function timestep(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh, equation = rhs
+
+ broadcast_volume!(cons2prim, equation, cache)
+ broadcast_volume!(speed, equation, cache)
+ max_v = get_variable(cache, :max_v)
+ vmax = dg1d.absolute_maximum(view(max_v,:))
+ vmax_limit = 0.99
+ if vmax > vmax_limit
+ @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
+ vmax = vmax_limit
+ end
+ smoothed_mu = get_variable(cache, :smoothed_mu)
+ mumax = dg1d.absolute_maximum(view(smoothed_mu, :))
+
+ @unpack dl, x, CFL, element = mesh
+ @unpack N = element
+
+ dt = CFL / (vmax * N / dl + mumax * N^3 / dl^2)
+ dt_limit = 1e-7
+ if dt < dt_limit
+ @warn "Limiting timestep due to dt being smaller than $dt_limit"
+ dt = dt_limit
+ end
+ return dt
+end
+
+
+#######################################################################
+# RHS computation #
+#######################################################################
+
+
+function rhs!(state_du, state_u, t, rhs::RHS, bdryconds::BoundaryConditions,
+ ldg_bdryconds::Maybe{BoundaryConditions}, av_bdryconds::Maybe{BoundaryConditions})
+
+ @unpack cache = rhs
+ wrap_dynamic_variables!(cache, state_u)
+ wrap_rhs_variables!(cache, state_du)
+ rhs!(rhs, t, rhs.hrsc, bdryconds, ldg_bdryconds, av_bdryconds)
+end
+
+
+function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
+
+ @unpack cache, mesh, equation, rsolver = rhs
+ @unpack D, S, tau = get_dynamic_variables(cache)
+ @unpack flx_D, flx_S, flx_tau = get_static_variables(cache)
+ @unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
+ @unpack lhs_numflx_D, rhs_numflx_D,
+ lhs_numflx_S, rhs_numflx_S,
+ lhs_numflx_tau, rhs_numflx_tau = get_bdry_variables(cache)
+
+ broadcast_volume!(cons2prim, equation, cache)
+ broadcast_volume!(flux, equation, cache)
+ broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+
+ compute_rhs_weak_form!(rhs_D, flx_D, lhs_numflx_D, rhs_numflx_D, mesh)
+ compute_rhs_weak_form!(rhs_S, flx_S, lhs_numflx_S, rhs_numflx_S, mesh)
+ compute_rhs_weak_form!(rhs_tau, flx_tau, lhs_numflx_tau, rhs_numflx_tau, mesh)
+
+ return
+end
+
+
+function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractReconstruction,
+ bdryconds, ldg_bdryconds, av_bdryconds)
+
+ @unpack cache, mesh, equation, rsolver = rhs
+ @unpack D, S, tau = get_dynamic_variables(cache)
+ @unpack flx_D, flx_S, flx_tau = get_static_variables(cache)
+ @unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
+ @unpack lhs_numflx_D, rhs_numflx_D,
+ lhs_numflx_S, rhs_numflx_S,
+ lhs_numflx_tau, rhs_numflx_tau = get_bdry_variables(cache)
+ @unpack flag = get_cell_variables(cache)
+
+ HRSC.reconstruct!(D, flag, hrsc, isperiodic=isperiodic(bdryconds), m=1e-10, limit_with_boundaries=false)
+ HRSC.reconstruct!(S, flag, hrsc, isperiodic=isperiodic(bdryconds), limit_with_boundaries=false)
+ HRSC.reconstruct!(tau, flag, hrsc, isperiodic=isperiodic(bdryconds), m=1e-10, limit_with_boundaries=false)
+
+ broadcast_volume!(cons2prim, equation, cache)
+ broadcast_volume!(flux, equation, cache)
+ broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
+
+ compute_rhs_weak_form!(rhs_D, flx_D, lhs_numflx_D, rhs_numflx_D, mesh)
+ compute_rhs_weak_form!(rhs_S, flx_S, lhs_numflx_S, rhs_numflx_S, mesh)
+ compute_rhs_weak_form!(rhs_tau, flx_tau, lhs_numflx_tau, rhs_numflx_tau, mesh)
+
+ return
+end
+
+
+function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
+ bdryconds, ldg_bdryconds, av_bdryconds)
+
+ @unpack cache, mesh, equation, rsolver,
+ ldg_rsolver, av_rsolver = rhs
+ @unpack D, S, tau = get_dynamic_variables(cache)
+ @unpack flx_D, flx_S, flx_tau,
+ ldg_D, ldg_S, ldg_tau,
+ flx_ldg_D, flx_ldg_S, flx_ldg_tau = get_static_variables(cache)
+ @unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
+ # TODO dispose of lhs/rhs numflx since we only need one value per interface
+ # Does a penalty flux violate this assumption?
+ @unpack lhs_numflx_D, rhs_numflx_D,
+ lhs_numflx_S, rhs_numflx_S,
+ lhs_numflx_tau, rhs_numflx_tau,
+ lhs_numflx_ldg_D, rhs_numflx_ldg_D,
+ lhs_numflx_ldg_S, rhs_numflx_ldg_S,
+ lhs_numflx_ldg_tau, rhs_numflx_ldg_tau = get_bdry_variables(cache)
+
+ ## solve auxiliary equation: q + ∂x u = 0
+ # improved syntax
+ # broadcast_volume!(ldg_flux, equation, grid)
+ # broadcast_faces!(mean_flux, ldg_rsolver, grid)
+ # broadcast_boundaryconditions!(mean_flux, ldg_bdryconds(t), grid)
+ broadcast_volume!(ldg_flux, equation, cache)
+ broadcast_faces!(mean_flux, ldg_rsolver, cache, mesh)
+ broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, t)
+ compute_rhs_weak_form!(ldg_D, flx_ldg_D, lhs_numflx_ldg_D, rhs_numflx_ldg_D, mesh)
+ compute_rhs_weak_form!(ldg_S, flx_ldg_S, lhs_numflx_ldg_S, rhs_numflx_ldg_S, mesh)
+ compute_rhs_weak_form!(ldg_tau, flx_ldg_tau, lhs_numflx_ldg_tau, rhs_numflx_ldg_tau, mesh)
+
+ ## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
+ broadcast_volume!(av_flux, equation, cache)
+ broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
+ broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, t)
+
+ broadcast_volume!(cons2prim, equation, cache)
+ broadcast_volume!(flux, equation, cache)
+ broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+
+ # TODO Remove doubled memeory and calls below by allowing to broadcast with accumalative operations
+ # e.g. broadcast_volume!([av_flux, flux], equation, cache)
+ # which only works if av_flux and flux use the same outputs
+ @. flx_D += flx_ldg_D
+ @. flx_S += flx_ldg_S
+ @. flx_tau += flx_ldg_tau
+ @. lhs_numflx_D += lhs_numflx_ldg_D
+ @. lhs_numflx_S += lhs_numflx_ldg_S
+ @. lhs_numflx_tau += lhs_numflx_ldg_tau
+ @. rhs_numflx_D += rhs_numflx_ldg_D
+ @. rhs_numflx_S += rhs_numflx_ldg_S
+ @. rhs_numflx_tau += rhs_numflx_ldg_tau
+
+ compute_rhs_weak_form!(rhs_D, flx_D, lhs_numflx_D, rhs_numflx_D, mesh)
+ compute_rhs_weak_form!(rhs_S, flx_S, lhs_numflx_S, rhs_numflx_S, mesh)
+ compute_rhs_weak_form!(rhs_tau, flx_tau, lhs_numflx_tau, rhs_numflx_tau, mesh)
+
+ return
+end
+
+
+
+### Desired new interface
+function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
+ bdryconds, ldg_bdryconds, av_bdryconds)
+
+ @unpack grid, equation, rsolver, ldg_rsolver, av_rsolver = rhs
+
+ @unpack D, S, tau = dynamic_variables(cache)
+ @unpack rhs_D, rhs_S, rhs_tau = rhs_variables(cache)
+ @unpack flx_D, flx_S, flx_tau,
+ ldg_D, ldg_S, ldg_tau, = static_variables(cache)
+ @unpack numflx_D, numflx_S, numflx_tau = bdry_variables(cache)
+
+ success = broadcast_volume!(cons2prim, equation, grid)
+
+ ## solve auxiliary equation: q + ∂x u = 0
+ # improved syntax
+ broadcast_volume!(ldg_flux, equation, grid)
+ broadcast_faces!(mean_flux, ldg_rsolver, grid)
+ broadcast_boundaryconditions!(mean_flux, ldg_bdryconds(t), grid)
+
+ compute_rhs_weak_form!(ldg_D, flx_D, numflx_D, mesh)
+ compute_rhs_weak_form!(ldg_S, flx_S, numflx_S, mesh)
+ compute_rhs_weak_form!(ldg_tau, flx_tau, numflx_tau, mesh)
+
+ ## compute rhs of regularized equation: ∂t u + ∂x (f + μ q) = 0
+ broadcast_volume!([av_flux, flux], equation, grid)
+ broadcast_faces!([mean_flux, lax_friedrich_flux], [av_rsolver, rsolver], grid)
+ broadcast_boundaryconditions!([mean_flux, lax_friedrich_flux],
+ [av_bdryconds(t), bdryconds(t)],
+ grid)
+
+ compute_rhs_weak_form!(rhs_D, flx_D, numflx_D, grid)
+ compute_rhs_weak_form!(rhs_S, flx_S, numflx_S, grid)
+ compute_rhs_weak_form!(rhs_tau, flx_tau, numflx_tau, grid)
+
+ return
+end
+
+### Other considerations
+# - @with_signature: What is the maximum tuple size one can pass around?
+# If reached should be swap over to SVec?
+# Or are normal vectors good enough?
+# - @with_signature: Think about @inline addition to macro.
diff --git a/src/refvector.jl b/src/refvector.jl
new file mode 100644
index 0000000000000000000000000000000000000000..11ce1e0fed0eb7229a6299e447204b59d6cf2633
--- /dev/null
+++ b/src/refvector.jl
@@ -0,0 +1,216 @@
+export RefVector, RefBlockVector, block, update!, eachblock
+
+
+#######################################################################
+# RefVector #
+#######################################################################
+
+
+"""
+ RefVector{T} <: DenseVector{T}
+
+A vector referencing a contiguous subrange of a `Vector{T}`.
+
+---
+
+ RefVector(v::Vector{T}[, offset, n])
+
+Construct a `RefVector` tracking `v` starting at `offset` and length `n`.
+Use `update!(refvec, vec[, offset, n])` to update the reference and/or subrange.
+
+---
+
+# Examples
+
+```julia
+julia> v = randn(100);
+
+julia> refvec = RefVector(v, 10, 50);
+
+julia> v[10] == refvec[1]
+true
+
+julia> length(refvec)
+50
+
+julia> update!(refvec, v, 20, 50);
+
+julia> v[10] == refvec[1]
+false
+```
+"""
+mutable struct RefVector{T} <: DenseVector{T}
+ ref::Base.RefValue{Vector{T}}
+ # so far, vec is actually optional, but it is included for testing to see whether
+ # the ref above can have any influence on performance
+ vec::Vector{T}
+ offset::Int64
+ n::Int64
+
+ function RefVector{T}(v::Vector{T}, offset::Integer, n::Integer) where T
+ checkrange(length(v), offset, n)
+ vec = unsafe_wrap(Vector{T}, pointer(v, offset), n)
+ return new{T}(Ref(v), vec, offset, n)
+ end
+end
+
+RefVector(v::Vector{T}, offset::Integer=1, n::Integer=length(v)) where T =
+ RefVector{T}(v, offset, n)
+RefVector(ref_v::Base.RefValue{Vector{T}}, offset::Integer=1, n::Integer=length(ref[])) where T =
+ RefVector(ref_v[], offset, n)
+
+# array interface
+@inline Base.size(rv::RefVector) = (rv.n,)
+@propagate_inbounds Base.getindex(rv::RefVector, i) = rv.vec[i]
+@propagate_inbounds Base.setindex!(rv::RefVector, val, i) = rv.vec[i] = val
+@inline Base.IndexStyle(A::RefVector) = IndexLinear()
+
+Base.unsafe_convert(::Type{Ptr{T}}, rv::RefVector{T}) where T = pointer(rv.ref[], rv.offset)
+# Base.unsafe_convert(::Type{Ptr{T}}, rv::RefVector{T}) where T = Base.unsafe_convert(Ptr{T}, rv.vec)
+
+# @propagate_inbounds Base.getindex(rv::RefVector, i::Integer) = rv.ref[][i+rv.offset-1]
+# @propagate_inbounds Base.setindex!(rv::RefVector, val, i::Integer) = rv.ref[][i+rv.offset-1] = val
+# Base.iterate(rv::RefVector) = iterate(rv.vec)
+# Base.iterate(rv::RefVector, state) = iterate(rv.vec, state)
+# Base.BroadcastStyle(::Type{RefVector{T}}) where T = Base.BroadcastStyle(Vector{T})
+# Base.similar(bc::Broadcasted{DestStyle}, ::Type{ElType})
+# Base.copyto!(dest, bc::Broadcasted{DestStyle}) =
+# function Base.show(io::IOBuffer, rv::RefVector)
+# println(io, "$(typeof(rv)):")
+# for v in rv
+# println(io, v)
+# end
+# return
+# end
+
+
+"""
+ update!(rv::RefVector{T}, v::Vector{T} [, offset, n]) where T
+
+Update the reference and/or subrange of `rv` to track `v`.
+"""
+function update!(rv::RefVector{T}, v::Vector{T}, offset=rv.offset, n=rv.n) where T
+ checkrange(length(v), offset, n)
+ rv.ref[] = v
+ rv.vec = unsafe_wrap(Vector{T}, pointer(v, offset), n)
+ rv.offset = offset
+ rv.n = n
+ return
+end
+
+
+function checkrange(N, offset, n)
+ if !(1 <= offset <= N-n+1) || !(0 <= n <= N)
+ error("Cannot reference vector of length $N starting from offset $offset with length $n")
+ end
+end
+
+
+#######################################################################
+# RefBlockVector #
+#######################################################################
+
+
+"""
+ RefBlockVector{T,T_Vector<:DenseVector{T}} <: DenseVector{T}
+
+A block vector referencing a `T_Vector<:DenseVector{T}`.
+
+---
+
+ RefBlockVector(v::Vector{T}, block_sizes[, offset])
+
+Construct a `RefBlockVector` tracking `v` by blocks of length `block_sizes`,
+optionally starting at `offset`.
+Use `update!(refvec, vec[, block_sizes, offset])` to update the reference and/or blocks.
+
+---
+
+# Examples
+
+```julia
+julia> v = randn(100);
+
+julia> refbvec = RefBlockVector(v, [10, 20, 30], 5)
+
+julia> v[15:25] == block(refbvec, 1)
+true
+
+julia> length(refbvec)
+60
+
+julia> update!(refvec, v, [10, 20], 1);
+
+julia> v[1:10] == block(refvec, 2)
+true
+```
+"""
+mutable struct RefBlockVector{T,T_Vector<:DenseVector{T}} <: DenseVector{T}
+ ref::Base.RefValue{T_Vector}
+ block_sizes::Vector{Int64}
+ blocks::Vector{RefVector{T}}
+ offset::Int64
+ n::Int64
+
+ function RefBlockVector{T,T_Vector}(v::T_Vector, block_sizes::Vector{Int64}, offset::Int64) where {T,T_Vector<:DenseVector{T}}
+ blocks = RefVector{T}[]
+ sizehint!(blocks, length(block_sizes))
+ refvec = new{T,T_Vector}(Ref(v), block_sizes, blocks)
+ update!(refvec, v, block_sizes, offset)
+ return refvec
+ end
+end
+
+
+RefBlockVector(v::T_Vector, block_sizes::Vector{Int64}, offset::Int64=1) where {T,T_Vector<:DenseVector{T}} =
+ RefBlockVector{T,T_Vector}(v, block_sizes, offset)
+
+
+"""
+ update!(refvec::RefBlockVector{T}, v::T_Vector, [block_sizes, offset]) where {T,T_Vector<:DenseVector{T}}
+
+Update the reference and/or blocks of `refvec` to track `v`.
+"""
+function update!(refvec::RefBlockVector{T}, v::T_Vector,
+ block_sizes::Vector{Int64}=refvec.block_sizes, offset::Int64=refvec.offset) where {T,T_Vector<:DenseVector{T}}
+
+ n = sum(block_sizes)
+ check_block_sizes(length(v), n, block_sizes, offset)
+ refvec.ref[] = v
+ refvec.block_sizes .= block_sizes
+ refvec.n = n
+ refvec.offset = offset
+ block_offset = offset
+ resize!(refvec.blocks, length(block_sizes))
+ for (i, bsz) in enumerate(block_sizes)
+ if isassigned(refvec.blocks, i)
+ update!(refvec.blocks[i], v, block_offset, bsz)
+ else
+ refvec.blocks[i] = RefVector{T}(v, block_offset, bsz)
+ end
+ block_offset += bsz
+ end
+
+ return
+end
+
+
+function check_block_sizes(N, n, block_sizes, offset)
+ if !(1 <= offset <= N-n+1) || !(0 <= n <= N)
+ error("Cannot reference vector of length $N starting from offset $offset with length $n")
+ end
+ if any(bsz -> bsz <= 0, block_sizes)
+ error("block_sizes must be positive")
+ end
+end
+
+
+# Abstract array interface
+@inline Base.size(vec::RefBlockVector) = (vec.n,)
+@propagate_inbounds Base.getindex(vec::RefBlockVector, i::Integer) = vec.ref[][i+vec.offset-1]
+@propagate_inbounds Base.setindex!(vec::RefBlockVector, val, i::Integer) = vec.ref[][i+vec.offset-1] = val
+@inline Base.IndexStyle(A::RefBlockVector) = IndexLinear()
+
+# utilities
+@propagate_inbounds block(vec::RefBlockVector, i) = vec.blocks[i]
+eachblock(vec::RefBlockVector) = (block(vec, i) for i = 1:length(vec.blocks))
diff --git a/src/visz_trees.jl b/src/visz_trees.jl
new file mode 100644
index 0000000000000000000000000000000000000000..2fadb6a3bd86460371933190a93592f71f46ab3a
--- /dev/null
+++ b/src/visz_trees.jl
@@ -0,0 +1,48 @@
+using GLMakie
+GLMakie.activate!()
+using dg1d
+
+
+function plot_tree1d()
+
+ t1d = Tree1d(5)
+ xs = [ Float64(c[1]) for c in t1d.coords ]
+ ys = zeros(length(xs))
+ is_bdry = dg1d.is_boundarycell.(t1d.cells)
+ bxs, bys = xs[is_bdry], ys[is_bdry]
+ fig = Figure()
+ scatterlines(fig[1,1], xs, ys)
+ scatter!(fig[1,1], bxs, bys, color=:red)
+
+ return fig
+end
+
+
+function plot_tree2d()
+
+ t2d = Tree2d(10,10)
+ xs_ys = [ t2d.coords[row][col] for row = 1:n_cells(t2d), col = 1:2 ]
+ xs, ys = xs_ys[:,1], xs_ys[:,2]
+ is_bdry = dg1d.is_boundarycell.(t2d.cells)
+ bxs, bys = xs[is_bdry], ys[is_bdry]
+ fig = Figure()
+ scatterlines(fig[1,1], xs, ys)
+ scatter!(fig[1,1], bxs, bys, color=:red)
+
+ return fig
+end
+
+
+function plot_tree3d()
+
+ t3d = Tree3d(10,10,10)
+ xs_ys_zs = [ t3d.coords[row][col] for row = 1:n_cells(t3d), col = 1:3 ]
+ xs, ys, zs = xs_ys_zs[:,1], xs_ys_zs[:,2], xs_ys_zs[:,3]
+ is_bdry = dg1d.is_boundarycell.(t3d.cells)
+ fig = Figure()
+ scatterlines(fig[1,1], xs, ys, zs, markersize=100)
+ bxs, bys, bzs = xs[is_bdry], ys[is_bdry], zs[is_bdry]
+ scatter!(fig[1,1], bxs, bys, bzs, color=:red, markersize=200)
+
+ return fig
+end
diff --git a/src/wip_dg_rhs.jl b/src/wip_dg_rhs.jl
new file mode 100644
index 0000000000000000000000000000000000000000..3a979ad524d104fd9603c2d0d2377399ee9be6b0
--- /dev/null
+++ b/src/wip_dg_rhs.jl
@@ -0,0 +1,188 @@
+export compute_rhs_weak_form!,
+ compute_rhs_strong_form!
+
+
+"""
+ compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh)
+ compute_rhs_weak_form!(rhs, f, s, nf_lhs, nf_rhs, mesh)
+"""
+function compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh::Mesh)
+ @unpack invjac, element = mesh
+ @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
+ @unpack K = mesh
+ shape = layout(mesh)
+ mat_rhs = reshape(view(rhs,:), shape)
+ mat_f = reshape(view(f,:), shape)
+ # mat_MDM = reshape(MDM, (Npts,Npts))
+ # mul! currently allocates, although, (I think) it shouldn't ...
+ # a = @allocated begin
+ mul!(mat_rhs, MDM, mat_f)
+ # end
+ # display(a)
+ @. mat_rhs -= (nf_rhs' * Ml_rhs - nf_lhs' * Ml_lhs)
+ mat_rhs .*= invjac
+ return
+end
+function compute_rhs_weak_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::Mesh)
+ compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh)
+ rhs .+= s
+ return
+end
+
+
+"""
+ compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh)
+ compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh)
+"""
+function compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh::Mesh)
+ @unpack invjac, element = mesh
+ @unpack Npts, invM, D, Ml_lhs, Ml_rhs = element
+ @unpack K = mesh
+ shape = layout(mesh)
+ mat_rhs = reshape(view(rhs,:), shape)
+ mat_f = reshape(view(f,:), shape)
+ # compute rhs = - D * f inplace to rhs
+ # benchmarking showed that splitting of the factor of -1 is faster
+ # because a temporary is not created then
+ mul!(mat_rhs, D, mat_f)
+ mat_rhs .*= -1
+ @. mat_rhs += (nf_rhs' * Ml_rhs - nf_lhs' * Ml_lhs)
+ mat_rhs .*= invjac
+ return
+end
+function compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::Mesh)
+ compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh)
+ rhs .+= s
+ return
+end
+
+
+function compute_rhs_weak_form!(rhs, sf1, sf2, nsf1, nsf2, grid::Grid2d)
+
+ for cell in grid.tree
+
+ basis = grid.bases[cell.index]
+ offset = grid.offsets[cell.index]
+
+ vol_rng = range_volume(basis, offset)
+ bdry_rng = range_boundary(basis, offset)
+
+ @views begin
+ vec_rhs = rhs[vol_idx]
+ vec_f1, vec_f2 = f1[vol_rng], f2[vol_rng]
+ # do we need four views here? I could imagine that two bigger ones suffice here to
+ # if we make use of the conservation property of fluxes
+ vec_nsf1, vec_nsf2 = nfs1[bdry_rng], nfs2[bdry_rng]
+
+ cell_rhs = reshape(vec_rhs, vol_sz)
+ cell_f1 = reshape(vec_f1, vol_sz)
+ cell_f2 = reshape(vec_f2, vol_sz)
+ # no reshaping needed because in 1D boundarys are lines
+ cell_nsf1 = vec_nsf1
+ cell_nsf2 = vec_nsf2
+ end
+
+ b1, b2 = basis
+ nx, ny = size(basis)
+
+ # add up derivatives
+ muladd!(cell_rhs1, b1.D, cell_f1, -1.0, false)
+ muladd!(transpose(cell_rhs1), cell_f2, b2.D, -1.0, 1.0)
+
+ # add up numerical fluxes
+ # xmin
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, 1, :)))
+ mat_rhs[I] += (sub_nf1max[i] * b1.Ml_rhs - sub_nf1min[i] * b1.Ml_lhs)
+ end
+ # xmax
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, nx, :)))
+ mat_rhs[I] += (sub_nf1max[i] * b1.Ml_rhs - sub_nf1min[i] * b1.Ml_lhs)
+ end
+ # ymin
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, :, 1)))
+ mat_rhs[I] += (sub_nf2max[i] * b2.Ml_rhs - sub_nf2min[i] * b2.Ml_lhs)
+ end
+ # ymax
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, :, ny)))
+ mat_rhs[I] += (sub_nf2max[i] * b2.Ml_rhs - sub_nf2min[i] * b2.Ml_lhs)
+ end
+
+ end
+ # apply jacobian
+ # TODO Benchmark whether this should be done here or inside the loop.
+ mat_rhs .*= invjac
+ return
+end
+function compute_rhs_weak_form!(rhs, s, sf1, sf2, nsf1, nsf2, grid::Grid2d)
+ compute_rhs_weak_form!(rhs, sf1, sf2, nsf1, nsf2, grid)
+ # an explicit source variable might not be needed, instead one could directly evaluate
+ # the source terms into the rhs when computing them
+ rhs .+= s
+ return
+end
+
+
+function compute_rhs_weak_form!(rhs, sf1, sf2, sf3, nsf1, nsf2, nsf3, grid::Grid3d)
+
+ for cell in grid.tree
+
+ basis = grid.bases[cell.index]
+ offset = grid.offsets[cell.index]
+
+ vol_rng = range_volume(basis, offset)
+ bdry_rng = range_boundary(basis, offset)
+
+ @views begin
+ vec_rhs = rhs[vol_idx]
+ vec_f1, vec_f2, vec_f3 = f1[vol_rng], f2[vol_rng], f3[vol_rng]
+ # do we need four views here? I could imagine that two bigger ones suffice here to
+ # if we make use of the conservation property of fluxes
+ vec_nsf1, vec_nsf2, vec_nsf3 = nfs1[bdry_rng], nfs2[bdry_rng], nsf3[bdry_rng]
+
+ cell_rhs = reshape(vec_rhs, vol_sz)
+ cell_f1 = reshape(vec_f1, vol_sz)
+ cell_f2 = reshape(vec_f2, vol_sz)
+ cell_f3 = reshape(vec_f3, vol_sz)
+ # no reshaping needed because in 1D boundarys are lines
+ cell_nsf1 = vec_nsf1
+ cell_nsf2 = vec_nsf2
+ end
+
+ b1, b2 = basis
+ nx, ny = size(basis)
+
+ # add up derivatives
+ muladd!(cell_rhs1, b1.D, cell_f1, -1.0, false)
+ muladd!(transpose(cell_rhs1), cell_f2, b2.D, -1.0, 1.0)
+
+ # add up numerical fluxes
+ # xmin
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, 1, :)))
+ mat_rhs[I] += (sub_nf1max[i] * b1.Ml_rhs - sub_nf1min[i] * b1.Ml_lhs)
+ end
+ # xmax
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, nx, :)))
+ mat_rhs[I] += (sub_nf1max[i] * b1.Ml_rhs - sub_nf1min[i] * b1.Ml_lhs)
+ end
+ # ymin
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, :, 1)))
+ mat_rhs[I] += (sub_nf2max[i] * b2.Ml_rhs - sub_nf2min[i] * b2.Ml_lhs)
+ end
+ # ymax
+ for (i, I) in enumerate(eachindex(@view(cell_rhs, :, ny)))
+ mat_rhs[I] += (sub_nf2max[i] * b2.Ml_rhs - sub_nf2min[i] * b2.Ml_lhs)
+ end
+
+ end
+ # apply jacobian
+ # TODO Benchmark whether this should be done here or inside the loop.
+ mat_rhs .*= invjac
+ return
+end
+function compute_rhs_weak_form!(rhs, s, sf1, sf2, sf3, nsf1, nsf2, nsf3, grid::Grid3d)
+ compute_rhs_weak_form!(rhs, sf1, sf2, sf3, nsf1, nsf2, nsf3, grid)
+ # an explicit source variable might not be needed, instead one could directly evaluate
+ # the source terms into the rhs when computing them
+ rhs .+= s
+ return
+end