text:
mutable struct CivecmI2Givens <: AbstractCivecm
    endogenous::Matrix{Float64}
    exogenous::Matrix{Float64}
    lags::Int64
    rank::Tuple{Int64,Int64}
    pars::Vector{Float64}
    llConvCrit::Float64
    maxiter::Int64
    Z0::Matrix{Float64}
    Z1::Matrix{Float64}
    Z2::Matrix{Float64}
    Z3::Matrix{Float64}
    R0::Matrix{Float64}
    R1::Matrix{Float64}
    R2::Matrix{Float64}
end

function civecmI2Givens(endogenous::Matrix{Float64}, exogenous::Matrix{Float64}, lags::Int64)
    mDX = diff(endogenous, dims = 1)
    mDDX = diff(diff(endogenous, dims = 1), dims = 1)
    mLDDX = lagmatrix(mDDX, 1:lags - 2)
    mDU = diff(exogenous, dims = 1)
    mDDU = diff(diff(exogenous, dims = 1), dims = 1)
    mLDDU = lagmatrix(mDDU, 0:lags - 2)

    Z0 = mDDX[lags - 1:end, :]
    Z1 = [mDX[lags - 1:end - 1, :] mDU[lags - 1:end - 1, :]]
    Z2 = [endogenous[lags:end - 1, :] exogenous[lags:end - 1, :]]
    Z3 = copy(mLDDX)
    if size(mLDDU, 1) > 0
        Z3 = [Z3 mLDDU]
    end
    if size(Z3, 2) > 0
        R0 = mreg(Z0, Z3)[2]
        R1 = mreg(Z1, Z3)[2]
        R2 = mreg(Z2, Z3)[2]
    else
        R0 = Z0
        R1 = Z1
        R2 = Z2
    end

    return CivecmI2Givens(endogenous, float64(exogenous), lags, (size(endogenous, 2), 0),
        Array{Float64}(npars(size(R1, 2), size(R0, 2), size(R0, 2), 0)), 1.0e-8, 50000, Z0, Z1, Z2, Z3, R0, R1, R2)
end
civecmI2Givens(endogenous::Matrix{Float64}, lags::Int64) = civecmI2Givens(endogenous, zeros(size(endogenous, 1), 0), lags)
civecmI2Givens(endogenous::Matrix{Float64}, exogenous::UnitRange, lags::Int64) = civecmI2Givens(endogenous, float64(reshape(exogenous, length(exogenous), 1)), lags)

function npars(p1::Integer,p::Integer,r::Integer,s::Integer)
    n = div(p*(p - 1),2) - div((p - r)*(p - r - 1),2) +
        div(p1*(p1 - 1),2) - div((p1 - r)*(p1 - r - 1),2) +
        r +
        p*r + r*(p1 - r) +
        div((p - r)*(p - r - 1),2) - div((p - r - s)*(p - r - s - 1),2) +
        div((p1 - r)*(p1 - r - 1),2) - div((p1 - r - s)*(p1 - r - s - 1),2) +
        s
    return n
end
npars(obj::CivecmI2Givens) = npars(size(obj.Z1, 2), size(obj.Z0, 2), obj.rank[1], obj.rank[2])

function setrank(obj::CivecmI2Givens, rank::Tuple{Int64, Int64})
    if sum(rank) > size(obj.endogenous, 2)
        error("Illegal choice of rank")
    else
        obj.rank = rank
    end
    return obj
end

residuals(obj::CivecmI2Givens) = obj.R0 - obj.R2*β(obj)*σ_α_β(obj)*α(obj)' - obj.R1*Γ(obj)'

function α(obj::CivecmI2Givens, full::Bool)
    p = size(obj.Z0, 2)
    count = 1
       ans = Matrix{Float64}(I, p, p)
    for i = 1:obj.rank[1]
        for j = i + 1:p
            Q = LinearAlgebra.Givens(i, j, cos(obj.pars[count]), sin(obj.pars[count]))
            mul!(Q, ans)
            count += 1
        end
    end
    return (full ? ans : ans[:,1:obj.rank[1]])
end
α(obj::CivecmI2Givens) = α(obj, false)

function β(obj::CivecmI2Givens, full::Bool)
    p = size(obj.Z0, 2)
    p1 = size(obj.Z1, 2)
    count = div(p*(p-1),2) - div((p-obj.rank[1])*(p-obj.rank[1]-1),2) + 1
    ans = Matrix{Float64}(I, p1, p1)
    for i = 1:obj.rank[1]
        for j = i + 1:p1
            Q = LinearAlgebra.Givens(i, j, cos(obj.pars[count]), sin(obj.pars[count]))
            mul!(Q, ans)
            count += 1
        end
    end
    return (full ? ans : ans[:,1:obj.rank[1]])
end
β(obj::CivecmI2Givens) = β(obj, false)

function σ_α_β(obj::CivecmI2Givens)
    p = size(obj.Z0, 2)
    p1 = size(obj.Z1, 2)
    count = div(p*(p-1),2)-div((p-obj.rank[1])*(p-obj.rank[1]-1),2) +
        div(p1*(p1-1),2)-div((p1-obj.rank[1])*(p1-obj.rank[1]-1),2) + 1
    return diagm(obj.pars[count:count + obj.rank[1] - 1])
end

function Γ(obj::CivecmI2Givens)
    p = size(obj.Z0, 2)
    p1 = size(obj.Z1, 2)
    count = div(p*(p-1),2)-div((p-obj.rank[1])*(p-obj.rank[1]-1),2) +
        div(p1*(p1-1),2)-div((p1-obj.rank[1])*(p1-obj.rank[1]-1),2) +
        obj.rank[1] + 1
    ans = Matrix{Float64}(undef, p, p1)
    ans[:,1:obj.rank[1]] = obj.pars[count:count + p*obj.rank[1] - 1]
    ans[1:obj.rank[1],obj.rank[1] + 1:end] = obj.pars[count + p*obj.rank[1]:count + p*obj.rank[1] + (p - obj.rank[1])*obj.rank[1] - 1]
    ans[obj.rank[1] + 1:end,obj.rank[1] + 1:end] = ξ(obj)*σ_ξ_η(obj)*ηpar(obj)'
    return α(obj, true)*ans*β(obj, true)'
end

function ξ(obj::CivecmI2Givens)
    p = size(obj.Z0, 2)
    p1 = size(obj.Z1, 2)
    count = div(p*(p-1),2)-div((p-obj.rank[1])*(p-obj.rank[1]-1),2) +
        div(p1*(p1-1),2)-div((p1-obj.rank[1])*(p1-obj.rank[1]-1),2) +
        obj.rank[1] +
        p*obj.rank[1] + (p - obj.rank[1])*obj.rank[1] + 1
    ans = Matrix{Float64}(I, p - obj.rank[1], obj.rank[2])
    for i = 1:obj.rank[2]
        for j = i + 1:p - obj.rank[1]
            Q = LinearAlgebra.Givens(i, j, cos(obj.pars[count]), sin(obj.pars[count]))
            mul!(Q, ans)
            count += 1
        end
    end
    return ans
end

function ηpar(obj::CivecmI2Givens)
    p = size(obj.Z0, 2)
    p1 = size(obj.Z1, 2)
    count = div(p*(p-1),2)-div((p-obj.rank[1])*(p-obj.rank[1]-1),2) +
        div(p1*(p1-1),2)-div((p1-obj.rank[1])*(p1-obj.rank[1]-1),2) +
        obj.rank[1] +
        p*obj.rank[1] + (p - obj.rank[1])*obj.rank[1] +
        div((p - obj.rank[1])*(p-obj.rank[1] - 1), 2) - div((p - sum(obj.rank))*(p - sum(obj.rank) - 1), 2) + 1
    ans = Matrix{Float64}(I, p1 - obj.rank[1], obj.rank[2])
    for i = 1:obj.rank[2]
        for j = i + 1:p1 - obj.rank[1]
            Q = LinearAlgebra.Givens(i, j, cos(obj.pars[count]), sin(obj.pars[count]))
            mul!(Q, ans)
            count += 1
        end
    end
    return ans
end

function σ_ξ_η(obj::CivecmI2Givens)
    p = size(obj.Z0, 2)
    p1 = size(obj.Z1, 2)
    count = div(p*(p-1),2)-div((p-obj.rank[1])*(p-obj.rank[1]-1),2) +
        div(p1*(p1-1),2)-div((p1-obj.rank[1])*(p1-obj.rank[1]-1),2) +
        obj.rank[1] +
        p*obj.rank[1] + (p - obj.rank[1])*obj.rank[1] +
        div((p - obj.rank[1])*(p-obj.rank[1] - 1), 2) - div((p - sum(obj.rank))*(p - sum(obj.rank) - 1), 2) +
        div((p1 - obj.rank[1])*(p1-obj.rank[1] - 1), 2) - div((p1 - sum(obj.rank))*(p1 - sum(obj.rank) - 1), 2) + 1
    return diagm(obj.pars[count:count + obj.rank[2] - 1])
end

function loglikelihood(obj::CivecmI2Givens, pars)
    obj.pars = pars
    ll = loglikelihood(obj)
    #println(σ_α_β(obj)[1,1], "\t", ll)
    return ll
end


text:
struct Bingo
    table::Matrix{String}
end

Bingo() = Bingo([cast() for _ in 1:5, _ in 1:5])

string(rand('a':'z', 10)...)
Base.setindex!(b::Bingo, v::String, inds...) = b.table[inds...] = v
Base.getindex(b::Bingo, inds...) = b.table[inds...]
Base.lastindex(b::Bingo) = lastindex(b.table)

save(b::Bingo) = save("bingo.pdf", b)

function save(filename, b::Bingo)
    name, ext = splitext(filename)
    Luxor.Drawing(595, 708, filename)
    Luxor.origin()

    tiles = Tiler(600, 650, 5, 5)
    translate(0, 30)
    background("#6a6703")
    for (pos, n) in tiles
        sethue("white")
        box(pos, tiles.tilewidth, tiles.tileheight, :fillpreserve)
        sethue("black")
        strokepath()
        fontface("Monaco")

        fsize = 18
        fontsize(fsize)
        lines = textlines(b[n], 120; rightgutter=3)
        if length(lines) > 5
            fsize = 12
            fontsize(fsize)
            lines = textlines(b[n], 120; rightgutter=3)
        end

        if n == 13
            sethue("red")
            fontsize(22)
            fontface("Helvetica Bold")
            # settext("<b>FREE</b>", pos - (0, 40), halign="center", markup=true)
            text("FREE", pos - (0, 40), halign=:center)
            fontface("Monaco")
            fontsize(fsize)
            textbox(lines, pos - (0, 32), leading=20, alignment=:center)
            sethue("black")
        else
            textbox(lines, pos - (0, 50), leading=20, alignment=:center)
        end
    end
    sethue("black")
    fontsize(35)
    text("P", Point(-245, -360), valign = :top)
    text("I", Point(-130, -360), valign = :top)
    text("Qu", Point(-20, -360), valign = :top)
    text("I", Point(105, -360), valign = :top)
    text("L", Point(220, -360), valign = :top)
    Luxor.finish()
end

function generate(;preview=true)
    b = Bingo()
    save(b)

    if preview
        Luxor.preview()
    end
end

function Base.show(io::IO, b::Bingo)
    println("Bingo sheet:")
    print(io, b[1])
    for each in b[2:end]
        println(io)
        print(io, each)
    end
end

ext: jl


text:
@doc raw"""
    postproc(model::Potts, param::Parameter, obs::MCObservableSet)

# Observables to be calculated
In the following, $m$ is total magnetization per site and $\epsilon$ is total energy per site.

- `"Binder Ratio"`
    - $R := \frac{\left \langle m^4 \right \rangle}{\left \langle m^2 \right\rangle^2}$
- `"Susceptibility"`
    - $\chi := \frac{N}{T}\left(\left\langle m^2\right\rangle\right)$
- `"Connected Susceptibility"`
    - $\frac{N}{T}\left(\left\langle m^2\right\rangle - \left\langle |m| \right\rangle^2\right)$
- `"Specific Heat"`
    - $\frac{N}{T^2}\left(\left\langle \epsilon^2\right\rangle - \left\langle \epsilon \right\rangle^2\right)$
"""
function postproc(model::Potts, param::Parameter, obs::MCObservableSet)
    nsites = numsites(model)
    T = convert(Float64, param["T"])
    beta = 1.0/T

    jk = jackknife(obs)
    jk["Binder Ratio"] = jk["Magnetization^4"] / (jk["Magnetization^2"]^2)
    jk["Susceptibility"] = (nsites*beta)*jk["Magnetization^2"]
    jk["Connected Susceptibility"] = (nsites*beta)*(jk["Magnetization^2"] - jk["|Magnetization|"]^2)
    jk["Specific Heat"] = (nsites*beta*beta)*(jk["Energy^2"] - jk["Energy"]^2)
    return jk
end


text:
module HC_wrapper

    using HomotopyContinuation: independent_normal
    using Base: get_uuid_name
    using Symbolics
    using HomotopyContinuation
    using DocStringExtensions
    using ..HarmonicBalance

    include("HC_wrapper/homotopy_interface.jl")
end


text:
function solve_s_nested_flexible(sj,inv_sigma,deltas,J,order,included,maxes,marketvars)
# Solve for s as p moves
if size(sj,1)>1
    sj = sj';
end
if marketvars!=nothing
    extra = size(marketvars,2);
else
    extra = 0;
end

s = sj;
bernO = convert.(Integer,order);
BERN = []
# @show J
for xj = 1:J
    if xj>1
        if xj<J
            perm_order = hcat(convert.(Integer,xj:J)', convert.(Integer,1:xj-1)');
            perm_order = convert.(Integer, perm_order)
        else
            perm_order = convert.(Integer,1:J)
            perm_order = perm_order .- 1
            perm_order[1] = J
        end
    else
        perm_order = convert.(Integer, 1:J)
    end
    s_perm = similar(s)
    for i = 1:J
        s_perm[:,i] = s[:,perm_order[i]];
    end
    # becomes design matrix
    BERN_xj = zeros(size(deltas,1),1)
# Market shares
for j = 1:J
     if included[xj,perm_order[j]]==1
        BERN_xj = [BERN_xj bern(s_perm[:,j], bernO[xj,1] )]
     end
end

# Market characteristics
if marketvars!=nothing
    for j = 1:size(marketvars,2)
        BERN_xj = [BERN_xj bern(s_perm[:,j], bern[xj,1])]
    end
end
BERN_xj = BERN_xj[:,2:end]
push!(BERN, BERN_xj)
end
out = []
for j = 1:J
    numVars = sum(included[j,:]) + extra;
    design = fullInteraction(BERN[j], convert(Integer,numVars), bernO[j,1])
    temp = deltas[1,j] .- design*inv_sigma[j]
    append!(out, temp)
end
return convert.(Float64, out)

end


text:
export
    AbstractLaneCoord,
    LaneCoord,
    LaneCoordIm

abstract AbstractLaneCoord <: AbstractTrackCoord
type LaneCoord <: AbstractLaneCoord
    trackid::Cint # id of the corresponding track
    s::Cdouble    # distance along the track
    t::Cdouble    # distance perpendicular to the track
    z::Cdouble    # height
    h::Cdouble    # heading [rad]
    p::Cdouble    # pitch [rad]
    r::Cdouble
    laneid::Cint
    offset::Cdouble

    LaneCoord() = new() # create one without setting any values
    LaneCoord(trackid::Integer, s::Real, t::Real, z::Real = 0.0, h::Real=0.0, p::Real=0.0, r::Real=0.0, laneid::Real=1, offset::Real=0.0 ) =
        new(trackid, s, t, z, h, p, r, laneid, offset)
end
immutable LaneCoordIm <: AbstractLaneCoord
    trackid::Cint # id of the corresponding track
    s::Cdouble    # distance along the track
    t::Cdouble    # distance perpendicular to the track
    z::Cdouble    # height
    h::Cdouble    # heading [rad]
    p::Cdouble    # pitch [rad]
    r::Cdouble
    laneid::Cint
    offset::Cdouble

    LaneCoordIm(trackid::Integer, s::Real, t::Real, z::Real = 0.0, h::Real=0.0, p::Real=0.0, r::Real=0.0, laneid::Real=1, offset::Real=0.0 ) =
        new(trackid, s, t, z, h, p, r, laneid, offset)
end
Base.convert(::Type{LaneCoordIm}, c::LaneCoord) = LaneCoordIm(c.trackid, c.s, c.t, c.z, c.h, c.p, c.r, c.laneid, c.offset)
Base.convert(::Type{LaneCoord}, c::LaneCoordIm) = LaneCoord(c.trackid, c.s, c.t, c.z, c.h, c.p, c.r, c.laneid, c.offset)

function ==(a::AbstractLaneCoord, b::AbstractLaneCoord)
    a.trackid == b.trackid &&
    isapprox(a.s, b.s) &&
    isapprox(a.t, b.t) &&
    isapprox(a.z, b.z) &&
    isapprox(a.h, b.h) &&
    isapprox(a.p, b.p) &&
    isapprox(a.r, b.r) &&
    a.laneid == b.laneid &&
    isapprox(a.offset, b.offset)
end

function show(io::IO, coord::AbstractLaneCoord)
    @printf(io, "(%d, %.16e, %.16e, %.16e, %.16e, %.16e, %.16e, %d, %.16e)",
                  coord.trackid, coord.s, coord.t, coord.z,
                  coord.h, coord.p, coord.r, coord.laneid, coord.offset)
end
function print(io::IO, coord::AbstractLaneCoord)
    @printf(io, "(%d, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %d, %.3f)",
                  coord.trackid, coord.s, coord.t, coord.z,
                  coord.h, coord.p, coord.r, coord.laneid, coord.offset)
end

copy(c::LaneCoord) = LaneCoord(c.trackid, c.s, c.t, c.z, c.h, c.p, c.r, c.laneid, c.offset)
copy(c::LaneCoordIm) = LaneCoordIm(c.trackid, c.s, c.t, c.z, c.h, c.p, c.r, c.laneid, c.offset)
deepcopy(coord::AbstractLaneCoord) = copy(coord)

function copy!(dest::LaneCoord, src::AbstractLaneCoord)
    dest.trackid = src.trackid
    dest.s = src.s
    dest.t = src.t
    dest.z = src.z
    dest.h = src.h
    dest.p = src.p
    dest.r = src.r
    dest.laneid = src.laneid
    dest.offset = src.offset
    dest
end

function plus!(a::LaneCoord, b::AbstractLaneCoord)
    (ccall( (:lanecoord_plusequal_lanecoord, LIB_ODRMGR), Void, (Ptr{Void},Ptr{Void}),
        pointer_from_objref(a), pointer_from_objref(b)))
    a
end
function plus!(a::LaneCoord, b::TrackCoord)
    (ccall( (:lanecoord_plusequal_trackcoord, LIB_ODRMGR), Void, (Ptr{Void},Ptr{Void}),
        pointer_from_objref(a), pointer_from_objref(b)))
    a
end

function init!(coord::LaneCoord)
    coord.trackid = 0
    coord.s = 0.0
    coord.t = 0.0
    coord.z = 0.0
    coord.h = 0.0
    coord.p = 0.0
    coord.r = 0.0
    coord.laneid = 0
    coord.offset = 0.0
    coord
end

max_stars_repo_name: sisl/OpenDriveManagerLite

max_issues_repo_name: sisl/OpenDriveManagerLite

max_forks_repo_name: sisl/OpenDriveManagerLite


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""IP address information for entries in the (plural) PodIPs field. Each entry includes:    IP: An IP address allocated to the pod. Routable at least within the cluster.

    IoK8sApiCoreV1PodIP(;
        ip=nothing,
    )

    - ip::String : ip is an IP address (IPv4 or IPv6) assigned to the pod
"""
mutable struct IoK8sApiCoreV1PodIP <: SwaggerModel
    ip::Any # spec type: Union{ Nothing, String } # spec name: ip

    function IoK8sApiCoreV1PodIP(;ip=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodIP, Symbol("ip"), ip)
        setfield!(o, Symbol("ip"), ip)
        o
    end
end # type IoK8sApiCoreV1PodIP

const _property_map_IoK8sApiCoreV1PodIP = Dict{Symbol,Symbol}(Symbol("ip")=>Symbol("ip"))
const _property_types_IoK8sApiCoreV1PodIP = Dict{Symbol,String}(Symbol("ip")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PodIP }) = collect(keys(_property_map_IoK8sApiCoreV1PodIP))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodIP }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodIP[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodIP }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodIP[property_name]

function check_required(o::IoK8sApiCoreV1PodIP)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodIP }, name::Symbol, val)
end


text:
using AbstractTrees

"Supertype for all nodes"
abstract type AbstractNode end

"""
    DataNode{D} <: AbstractNode

Holds literal data of type D
"""
struct DataNode{D} <: AbstractNode
    data::D
end

Base.:(==)(n1::T, n2::T) where {T<:DataNode} = n1.data == n2.data

AbstractTrees.printnode(io::IO, node::DataNode) = print(io, node.data)

"""
    Node{T} <: Abstract Node

Holds children as vector. Type parameter `T` is a
symbol to represent the node type.
"""
mutable struct Node{T} <: AbstractNode
    children::Vector{AbstractNode}
    Node(type) = new{type}(AbstractNode[])
end

Base.:(==)(n1::T, n2::T) where {T<:Node} = n1.children == n2.children

type(::Node{T}) where {T} = T

AbstractTrees.children(node::Node) = node.children
AbstractTrees.printnode(io::IO, node::Node) = print(io, type(node))

function Base.match(stream::TokenStream, type::Symbol;  kwargs...)
    children = AbstractNode[]
    match!(children, stream, type; kwargs...)
end

function match!(children, stream::TokenStream, type::Symbol; needsmatch=true)
    candidate = Node(type)
    pos = position(stream)
    try
        construct(candidate, stream)
        push!(children, candidate)
        return candidate
    catch err
        if err isa ParsingError && !needsmatch
            setposition!(stream, pos) # reset! stream position
        else
            rethrow()
        end
    end
    return nothing
end

function match_many(args...; kwargs...)
    children = AbstractNode[]
    match_many!(children, args...; kwargs...)
    return children
end

function match_many!(children, stream, type; needsmatch=true, delimiter=nothing)
    counter = 0
    while true
        # try to match type, if its the first one check needsmatch
        found = match!(children, stream, type;
                       needsmatch = (needsmatch && iszero(counter)))
        found === nothing && break # break if no match
        counter += 1

        if delimiter !== nothing
            delim = match_token(stream, delimiter, needsmatch=false)
            # XXX is this right? break scan if no delim found
            delim === nothing && break
        end
    end
    return counter
end

function match_token(stream::TokenStream, type::Symbol, lexme=nothing; needsmatch=true)
    t = peek(stream)

    if t!==nothing && type === t.type && (lexme === nothing  || string(lexme) == t.lexme)
        return next!(stream)
    else
        if needsmatch
            throw(ParsingError("Wrong token! got '$t' but requested '$(Token(type, lexme))'", stream))
        end
        return nothing
    end
end

function construct(n::Node{:bool_literal}, stream)
    t = match_token(stream, :bool_literal)
    d = DataNode(parse(Bool, t.lexme))
    n.children = [d]
end

function construct(n::Node{:int_literal}, stream)
    t = match_token(stream, :int_literal)
    d = DataNode(parse(Int, t.lexme))
    n.children = [d]
end

function construct(n::Node{:float_literal}, stream)
    t = match_token(stream, :float_literal)
    d = DataNode(parse(Float64, t.lexme))
    n.children = [d]
end

function construct(n::Node{:string_literal}, stream)
    t = match_token(stream, :string_literal)
    d = DataNode(t.lexme)
    n.children = [d]
end

function construct(n::Node{:identifier}, stream)
    t = match_token(stream, :identifier)
    d = DataNode(t.lexme)
    n.children = [d]
end

function construct(n::Node{:model}, stream)
    match_many!(n.children, stream, :predicate_item, needsmatch=false)
    match_many!(n.children, stream, :par_decl_item, needsmatch=false)
    match_many!(n.children, stream, :var_decl_item, needsmatch=false)
    match_many!(n.children, stream, :constraint_item, needsmatch=false)
    match!(n.children, stream, :solve_item)
end

function construct(n::Node{:fullline}, stream)
    if match!(n.children, stream, :constraint_item, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :predicate_item, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :par_decl_item, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :var_decl_item, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :solve_item) !== nothing
        return
    end
    throw(ParsingError("Could not construct :fullline", stream))
end

function construct(n::Node{:predicate_item}, stream)
    match_token(stream, :keyword, "predicate")
    match!(n.children, stream, :identifier)
    match_token(stream, :parentheses_l)

    while true
        match!(n.children, stream, :pred_param_type)
        match_token(stream, :colon)
        match!(n.children, stream, :identifier)

        comma = match_token(stream, :comma, needsmatch=false)
        if comma === nothing
            match_token(stream, :parentheses_r)
            break
        end
        closing = match_token(stream, :parentheses_r, needsmatch=false)
        if closing !== nothing
            # found closing delim
            break
        end
    end
    match_token(stream, :semicolon)
end

function construct(n::Node{:pred_param_type}, stream)
    if match!(n.children, stream, :basic_pred_param_type, needsmatch=false) !== nothing
        return
    elseif match_token(stream, :keyword, "array", needsmatch=false) !== nothing
        match_token(stream, :bracket_l)
        match!(n.children, stream , :pred_index_set)
        match_token(stream, :bracket_r)
        match_token(stream, :keyword, "of")
        match!(n.children, stream, :basic_pred_param_type)
        return
    end
    throw(ParsingError("Could not construct :pred_param_type", stream))
end

function construct(n::Node{:pred_index_set}, stream)
    t = match_token(stream, :basic_par_type, "int")
    if t  !== nothing
        d = DataNode(t.lexme)
        n.children = [d]
        return
    elseif match!(n.children, stream, :index_set, needsmatch=false) !==nothing
        return
    end
    throw(ParsingError("Could not construct :pred_index_set", stream))
end

function construct(n::Node{:index_set}, stream)
    match_token(stream, :int_literal, "1")
    match_token(stream, :dotdot)
    match!(n.children, stream, :int_literal)
end

function construct(n::Node{:basic_pred_param_type}, stream)
    # try to match basic_par_type
    if match!(n.children, stream, :basic_par_type, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :basic_var_type, needsmatch=false) !== nothing
        return
    end

    # jump over "set of" if given
    if match_token(stream, :keyword, "set", needsmatch=false) !== nothing
        match_token(stream, :keyword, "of")
        # but catch "set of int" type that lexer does not support (currently)
        # XXX: more elegant solution would be to catch that in lexer,
        # but with the spaces this is somewhat tricky
        if match_token(stream, :basic_par_type, "int", needsmatch=false) !== nothing
            d = DataNode("set of int")
            push!(n.children, d)
            return
        end
    end

    if match!(n.children, stream, :set_literal, needsmatch=false) !== nothing
        # XXX: this will also match {Float, FLoat} which is not in the grammar
        return
    end

    throw(ParsingError("Could not construct :basic_pred_param_type", stream))
end

function construct(n::Node{:basic_par_type}, stream)
    t = match_token(stream, :basic_par_type)
    d = DataNode(t.lexme)
    n.children = [d]
end

function construct(n::Node{:basic_var_type}, stream)
    match_token(stream, :keyword, "var")

    if match!(n.children, stream, :basic_par_type, needsmatch=false) !== nothing
        return
    end

    # jump over "set of" if given
    if match_token(stream, :keyword, "set", needsmatch=false) !== nothing
        match_token(stream, :keyword, "of")
    end

    if match!(n.children, stream, :set_literal, needsmatch=false) !== nothing
        # XXX: this will also match {Float, FLoat} which is not in the grammar
        return
    end

    throw(ParsingError("Could not construct :basic_par_type", stream))
end

function construct(n::Node{:set_literal}, stream)
    if match_token(stream, :brace_l, needsmatch=false) !== nothing
        if peek(stream).type == :int_literal
            match_many!(n.children, stream, :int_literal, delimiter=:comma)
        elseif peek(stream).type == :float_literal
            match_many!(n.children, stream, :float_literal, delimiter=:comma)
        else
            throw(ParsingError("Could not parse set_literal!", stream))
        end
        match_token(stream, :brace_r)
        return
    elseif match!(n.children, stream, :range, needsmatch=true) !== nothing
        return
    end
    throw(ParsingError("Could not construct :set_literal", stream))
end

function construct(n::Node{:range}, stream)
    if match!(n.children, stream, :int_literal, needsmatch=false) !== nothing
        match_token(stream, :dotdot)
        match!(n.children, stream, :int_literal)
        return
    elseif match!(n.children, stream, :float_literal, needsmatch=false) !== nothing
        match_token(stream, :dotdot)
        match!(n.children, stream, :float_literal)
        return
    end
    throw(ParsingError("Could not construct :range", stream))
end

function construct(n::Node{:par_decl_item}, stream)
    match!(n.children, stream, :par_type)
    match_token(stream, :colon)
    match!(n.children, stream, :identifier)
    match_token(stream, :equalsign)
    match!(n.children, stream, :par_expr)
    match_token(stream, :semicolon)
end

function construct(n::Node{:par_type}, stream)
    if match!(n.children, stream, :basic_par_type, needsmatch=false) !== nothing
        return
    elseif match_token(stream, :keyword, "array", needsmatch=false) !== nothing
        match_token(stream, :bracket_l)
        match!(n.children, stream, :index_set)
        match_token(stream, :bracket_r)
        match_token(stream, :keyword, "of")
        match!(n.children, stream, :basic_par_type)
        return
    end
    throw(ParsingError("Could not construct :par_type", stream))
end

function construct(n::Node{:par_expr}, stream)
    if match!(n.children, stream, :basic_literal_expr, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :par_array_literal, needsmatch=false) !== nothing
        return
    end
    throw(ParsingError("Could not construct :par_expr", stream))
end

function construct(n::Node{:basic_literal_expr}, stream)
    if match!(n.children, stream, :set_literal, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :bool_literal, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :int_literal, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :float_literal, needsmatch=false) !== nothing
        return
    end
    throw(ParsingError("Could not construct :basic_literal_expr", stream))
end

function construct(n::Node{:par_array_literal}, stream)
    match_token(stream, :bracket_l)
    match_many!(n.children, stream, :basic_literal_expr, delimiter=:comma)
    match_token(stream, :bracket_r)
end

function construct(n::Node{:var_decl_item}, stream)
    if match!(n.children, stream, :basic_var_type, needsmatch=false) !== nothing
        match_token(stream, :colon)
        match!(n.children, stream, :identifier)
        # well the type annotions is required but might be empty
        match!(n.children, stream, :annotations, needsmatch=false)
        if match_token(stream, :equalsign, needsmatch=false) !== nothing
            match!(n.children, stream, :par_expr)
        end
        match_token(stream, :semicolon)
        return
    elseif match!(n.children, stream, :array_var_type, needsmatch=false) !== nothing
        match_token(stream, :colon)
        match!(n.children, stream, :identifier)
        # well the type annotions is required but might be empty
        match!(n.children, stream, :annotations, needsmatch=false)
        match_token(stream, :equalsign)
        match!(n.children, stream, :array_literal)
        match_token(stream, :semicolon)
        return
    end
    throw(ParsingError("Could not construct :var_decl_item", stream))
end

function construct(n::Node{:array_var_type}, stream)
    match_token(stream, :keyword, "array")
    match_token(stream, :bracket_l)
    match!(n.children, stream, :index_set)
    match_token(stream, :bracket_r)
    match_token(stream, :keyword, "of")
    match!(n.children, stream, :basic_var_type)
end

function construct(n::Node{:array_literal}, stream)
    # needs to start with [
    match_token(stream, :bracket_l)
    match_many!(n.children, stream, :basic_expr, needsmatch=false, delimiter=:comma)
    match_token(stream, :bracket_r)
end

function construct(n::Node{:annotations}, stream)
    # needs to start with doublecolon
    match_token(stream, :doublecolon, needsmatch=true)
    match_many!(n.children, stream, :annotation; delimiter=:doublecolon, needsmatch=true)
end

function construct(n::Node{:annotation}, stream)
    # needs to start with identifier
    match!(n.children, stream, :identifier)
    if match_token(stream, :parentheses_l, needsmatch=false) !== nothing
        match_many!(n.children, stream, :ann_expr, delimiter=:comma)
        match_token(stream, :parentheses_r)
    end
end

function construct(n::Node{:ann_expr}, stream)
    if match!(n.children, stream, :basic_ann_expr, needsmatch=false) !== nothing
        return
    elseif match_token(stream, :bracket_l, needsmatch=false) !== nothing
        match_many!(n.children, stream, :basic_ann_expr, delimiter=:comma)
        match_token(stream, :bracket_r)
        return
    end
    throw(ParsingError("Could not construct :ann_expr", stream))
end

function construct(n::Node{:basic_ann_expr}, stream)
    if match!(n.children, stream, :annotation, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :basic_literal_expr, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :string_literal, needsmatch=false) !== nothing
        return
    end
    throw(ParsingError("Could not construct :basic_ann_expr", stream))
end

function construct(n::Node{:expr}, stream)
    if match!(n.children, stream, :basic_expr, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :array_literal, needsmatch=false) !== nothing
        return
    end
    throw(ParsingError("Could not construct :expr", stream))
end

function construct(n::Node{:basic_expr}, stream)
    if match!(n.children, stream, :basic_literal_expr, needsmatch=false) !== nothing
        return
    elseif match!(n.children, stream, :identifier, needsmatch=false) !== nothing
        return
    end
    throw(ParsingError("Could not construct :basic_expr", stream))
end

function construct(n::Node{:constraint_item}, stream)
    match_token(stream, :keyword, "constraint")
    match!(n.children, stream, :identifier)
    match_token(stream, :parentheses_l)
    match_many!(n.children, stream, :expr, delimiter=:comma)
    match_token(stream, :parentheses_r)
    match!(n.children, stream, :annotations, needsmatch=false) #technically nm=true
    match_token(stream, :semicolon)
end

function construct(n::Node{:solve_item}, stream)
    match_token(stream, :keyword, "solve")
    match!(n.children, stream, :annotations, needsmatch=false) #technically nm=true
    if match_token(stream, :keyword, "satisfy", needsmatch=false) !== nothing
        # nothing else to match
    elseif match_token(stream, :keyword, "minimize", needsmatch=false) !== nothing
        match!(n.children, stream, :basic_expr)
    elseif match_token(stream, :keyword, "maximize", needsmatch=false) !== nothing
        match!(n.children, stream, :basic_expr)
    else
        throw(ParsingError("Could not construct :solve_item", stream))
    end
    match_token(stream, :semicolon)
end


text:
function rhs!(du, u, t,
              mesh::CurvedMesh{2}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
    @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # solve the auxilary system for viscous flows
    if typeof(equations) == CompressibleEulerEquations2D{Float64} && equations.viscous
        #  q is an array with one dimension more than u: q = [∂u/∂x] [∂u/∂y]
        @unpack q1 = cache
        @unpack q2 = cache
        q1 .= zero(eltype(q1))
        q2 .= zero(eltype(q2))
        @timeit_debug timer() "auxiliary system" calc_nabla_u!(q1,
        q2, u, t, mesh, Trixi.AuxiliaryEquation(),
        boundary_conditions, # hier noch unstetige boundaries schwierig
        dg, cache)
        #println(maximum(abs.(q1)))
        # q1 .= zero(eltype(q1))
        # q2 .= zero(eltype(q2))
    end

  # Calculate volume integral
    @timeit_debug timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Calculate interface fluxes
    @timeit_debug timer() "interface flux" calc_interface_flux!(
    cache, u, mesh, equations, dg)

  # Calculate boundary fluxes
    @timeit_debug timer() "boundary flux" calc_boundary_flux!(
    cache, u, t, boundary_conditions, mesh, equations, dg)

  # Calculate surface integrals
    @timeit_debug timer() "surface integral" calc_surface_integral!(
    du, mesh, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
    @timeit_debug timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
    @timeit_debug timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

    return nothing
end

# calculation of the auxilary system q = ∇u for viscous flows
function calc_nabla_u!(q1, q2, u, t,
    mesh::CurvedMesh{2}, equations,
    boundary_conditions,
    dg::DG, cache)
    # use the standard weak form DGSEM to calculate ∇u
    # Calculate volume integral
    calc_volume_integral_auxiliary!(q1, q2, u,
        mesh::CurvedMesh{2}, equations::AuxiliaryEquation,
        dg::DG, cache)

    # calculate surface integral
    calc_surface_integral_auxiliary!(q1, q2, u, t,
        mesh::CurvedMesh{2}, equations::AuxiliaryEquation,
        dg::DG, cache,
        boundary_conditions)

  # Apply Jacobian from mapping to reference element
    apply_jacobian!(
        q1, mesh, equations, dg, cache)
    apply_jacobian!(
        q2, mesh, equations, dg, cache)

    return nothing
end

function calc_surface_integral_auxiliary!(q1, q2, u, t,
                                          mesh::CurvedMesh{2}, equations::AuxiliaryEquation,
                                          dg::DG, cache,
                                          boundary_conditions)
  # values are stored in surface_flux_values and calculated with u; then copied to q
    @unpack elements = cache
  # calc surface integral for q1
    @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2)"

    # Interfaces for u_x (1 dispatch) and in x-direction (`orientation` = 1)
        calc_interface_flux_auxiliary!(1, elements.surface_flux_values,
                         elements.left_neighbors[1, element],
                         element, 1, u, mesh, equations, dg, cache)

    # Interfaces Interfaces for u_x (1 dispatch) and in y-direction (`orientation` = 2)
        calc_interface_flux_auxiliary!(1, elements.surface_flux_values,
                         elements.left_neighbors[2, element],
                         element, 2, u, mesh, equations, dg, cache)
    end
  # boundary flux with nabla = 1 dispatch for u_x
    calc_boundary_flux_auxiliary!(cache, u, t, boundary_conditions, mesh, equations, dg, 1)
  # add the surface integral to q1
    calc_surface_integral!(q1, mesh, equations, dg, cache)

  # calc surface integral for q2. The previously calculated surface fluxes are overwritten.
    @threaded for element in eachelement(dg, cache)
        # Interfaces in negative directions
        # Faster version of "for orientation in (1, 2)"

        # Interfaces Interfaces for u_y (2 dispatch) and in x-direction (`orientation` = 1)
        calc_interface_flux_auxiliary!(2, elements.surface_flux_values,
                             elements.left_neighbors[1, element],
                             element, 1, u, mesh, equations, dg, cache)

        # Interfaces Interfaces for u_y (2 dispatch) and in y-direction (`orientation` = 2)
        calc_interface_flux_auxiliary!(2, elements.surface_flux_values,
                             elements.left_neighbors[2, element],
                             element, 2, u, mesh, equations, dg, cache)
    end
    # boundary flux with nabla = 2 dispatch for u_y
    calc_boundary_flux_auxiliary!(cache, u, t, boundary_conditions, mesh, equations, dg, 2)
    # add the surface integral to q2
    calc_surface_integral!(q2, mesh, equations, dg, cache)

    return nothing
end

function calc_volume_integral_auxiliary!(q1, q2, u,
                                        mesh::CurvedMesh{2}, equations::AuxiliaryEquation,
                                        dg::DGSEM, cache)
    @unpack derivative_dhat = dg.basis
    @unpack contravariant_vectors = cache.elements

    @threaded for element in eachelement(dg, cache)
        for j in eachnode(dg), i in eachnode(dg)
            # flux(u) = - u !
            u_node = - get_node_vars(u, equations, dg, i, j, element)

      # Compute the contravariant by taking the scalar product of the
      # first contravariant vector Ja^1 and u
            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
            contravariant_u_x1 = Ja11 * u_node
            contravariant_u_y1 = Ja12 * u_node
            for ii in eachnode(dg)
                integral_contribution1 = derivative_dhat[ii, i] * contravariant_u_x1
                integral_contribution2 = derivative_dhat[ii, i] * contravariant_u_y1
                add_to_node_vars!(q1, integral_contribution1, equations, dg, ii, j, element)
                add_to_node_vars!(q2, integral_contribution2, equations, dg, ii, j, element)
            end

      # Compute the contravariant by taking the scalar product of the
      # second contravariant vector Ja^2 and u
            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
            contravariant_u_x2 = Ja21 * u_node
            contravariant_u_y2 = Ja22 * u_node

            for jj in eachnode(dg)
                integral_contribution1 = derivative_dhat[jj, j] * contravariant_u_x2
                integral_contribution2 = derivative_dhat[jj, j] * contravariant_u_y2
                add_to_node_vars!(q1, integral_contribution1, equations, dg, i, jj, element)
                add_to_node_vars!(q2, integral_contribution2, equations, dg, i, jj, element)
            end
        end
    end
    return nothing
end

function calc_volume_integral!(du, u,
                               mesh::CurvedMesh{2},
                               nonconservative_terms::Val{false}, equations,
                               volume_integral::VolumeIntegralWeakForm,
                               dg::DGSEM, cache)
    @unpack derivative_dhat = dg.basis
    @unpack contravariant_vectors = cache.elements

    @threaded for element in eachelement(dg, cache)
        for j in eachnode(dg), i in eachnode(dg)
            u_node = get_node_vars(u, equations, dg, i, j, element)

            flux1 = flux(u_node, 1, equations)
            flux2 = flux(u_node, 2, equations)

      # Compute the contravariant flux by taking the scalar product of the
      # first contravariant vector Ja^1 and the flux vector
            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
            contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2

            for ii in eachnode(dg)
                integral_contribution = derivative_dhat[ii, i] * contravariant_flux1
                add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
            end

      # Compute the contravariant flux by taking the scalar product of the
      # second contravariant vector Ja^2 and the flux vector
            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
            contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2

            for jj in eachnode(dg)
                integral_contribution = derivative_dhat[jj, j] * contravariant_flux2
                add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
            end
        end
    end
    return nothing
end

##################BLZ: Volume Integral Pseudo Strong Form ##############
function calc_volume_integral!(du, u,
    mesh::CurvedMesh{2},
    nonconservative_terms::Val{false}, equations,
    volume_integral::VolumeIntegralPseudoStrongForm,
    dg::DGSEM, cache)
    @unpack derivative_matrix = dg.basis
    @unpack contravariant_vectors = cache.elements
    @unpack weights = dg.basis
  # Estimate the Pseudo strong Form D_pseudo = - M^-1 B + D
    D_pseudo = derivative_matrix
    D_pseudo[1, 1] += 1 / weights[1]
    D_pseudo[end,end] -= 1 / weights[end]

    @threaded for element in eachelement(dg, cache)
        for j in eachnode(dg), i in eachnode(dg)
            u_node = get_node_vars(u, equations, dg, i, j, element)

            flux1 = flux(u_node, 1, equations)
            flux2 = flux(u_node, 2, equations)

  # Compute the contravariant flux by taking the scalar product of the
  # first contravariant vector Ja^1 and the flux vector
            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
            contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2

            for ii in eachnode(dg)
                integral_contribution = D_pseudo[ii, i] * contravariant_flux1
                add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
            end

    # Compute the contravariant flux by taking the scalar product of the
    # second contravariant vector Ja^2 and the flux vector
            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
            contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2

            for jj in eachnode(dg)
                integral_contribution = D_pseudo[jj, j] * contravariant_flux2
                add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
            end
        end
    end
    return nothing
end

##################BLZ: Volume Integral FLuxDifferencing based on PseudoStrong ##############
function calc_volume_integral!(du, u,
    mesh::CurvedMesh{2},
    nonconservative_terms::Val{false}, equations,
    volume_integral::VolumeIntegralFluxDifferencing,
    dg::DGSEM, cache)

    if typeof(equations) == CompressibleEulerEquations2D{Float64} && equations.viscous
        calc_volume_integral_viscous!(du, u,
        mesh::CurvedMesh{2},
        nonconservative_terms::Val{false}, equations,
        volume_integral::VolumeIntegralFluxDifferencing,
        dg::DGSEM, cache)
    else
    # separate and parallel computation for each element
        @threaded for element in eachelement(dg, cache)
            @unpack contravariant_vectors = cache.elements
            @unpack inverse_weights = dg.basis
            @unpack derivative_matrix = dg.basis
            D_pseudo = derivative_matrix
            D_pseudo[1, 1] += inverse_weights[1]
            D_pseudo[end,end] -= inverse_weights[end]

    # Volume Integral for one Element, loop over all nodes
            for j in eachnode(dg), i in eachnode(dg)
                u_node = get_node_vars(u, equations, dg, i, j, element)
                flux1 = flux(u_node, 1, equations)
                flux2 = flux(u_node, 2, equations)

    # x direction
                Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)

    # two point flux of the same node outside the loop. And Pseudo Strong Part for boundary nodes
                if i == 1 || i == nnodes(dg)
                # Here we can also use the phys. flux.
                    contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
                    integral_contribution = (2 * D_pseudo[i, i] + S[i, i]) *
                                        contravariant_flux1
                    add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
                else
                    contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
                    integral_contribution = 2 * D_pseudo[i, i] * contravariant_flux1
                    add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
                end

    # two point flux symmetry: estimate the volume flux only one time
                for ii in (i + 1):nnodes(dg)
                # arihtmetic mean for two point volume flux metrics
                    Ja11_, Ja12_ = get_contravariant_vector(1, contravariant_vectors, ii, j, element)
                    u_node_ii = get_node_vars(u, equations, dg, ii, j, element)
                    flux1 = volume_integral.volume_flux(u_node, u_node_ii, 1, equations)
                    flux2 = volume_integral.volume_flux(u_node, u_node_ii, 2, equations)
                    contravariant_flux1 = 0.5 * (Ja11 + Ja11_) * flux1 + 0.5 * (Ja12 + Ja12_) * flux2
                    integral_contribution = 2 * D_pseudo[i, ii] * contravariant_flux1
                    add_to_node_vars!(du, integral_contribution, equations, dg, i,  j, element)
                    integral_contribution = 2 * D_pseudo[ii, i] * contravariant_flux1
                    add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
                end
    # y direction
                flux1 = flux(u_node, 1, equations)
                flux2 = flux(u_node, 2, equations)

                Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)

    # two point flux of the same node outside the loop. And Pseudo Strong Part for boundary nodes
                if j == 1 || j == nnodes(dg)
                # Here we can also use the phys. flux.
                    contra_flux = Ja22 * flux2 + Ja21 * flux1
                    integral_contribution = (2 * D_pseudo[j, j] + S[j, j]) *
                                         contra_flux
                    add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
                else
                    contravariant_flux2 = Ja22 * flux2 + Ja21 * flux1
                    integral_contribution = 2 * D_pseudo[j, j] * contravariant_flux2
                    add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
                end
    # two point flux symmetry: estimate the volume flux only one time
                for jj in (j + 1):nnodes(dg)
                # arithmetic mean for two point volume flux metrics
                    Ja21_, Ja22_ = get_contravariant_vector(2, contravariant_vectors, i, jj, element)
                    u_node_jj = get_node_vars(u, equations, dg, i, jj, element)
                    flux1 = volume_integral.volume_flux(u_node, u_node_jj, 1, equations)
                    flux2 = volume_integral.volume_flux(u_node, u_node_jj, 2, equations)
                    contravariant_flux2 = 0.5 * (Ja22 + Ja22_) * flux2 + 0.5 * (Ja21 + Ja21_) * flux1
                    integral_contribution =  2 * D_pseudo[j, jj] * contravariant_flux2
                    add_to_node_vars!(du, integral_contribution, equations, dg, i, j,  element)
                    integral_contribution =  2 * D_pseudo[jj, j] * contravariant_flux2
                    add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
                end
            end
        end
    end
end

function calc_volume_integral_viscous!(du, u,
    mesh::CurvedMesh{2},
    nonconservative_terms::Val{false}, equations,
    volume_integral::VolumeIntegralFluxDifferencing,
    dg::DGSEM, cache)

    @unpack q1 = cache
    @unpack q2 = cache

    # separate and parallel computation for each element
    @threaded for element in eachelement(dg, cache)
        @unpack contravariant_vectors = cache.elements
        @unpack inverse_weights = dg.basis
        @unpack derivative_matrix = dg.basis
        D_pseudo = derivative_matrix
        D_pseudo[1, 1] += inverse_weights[1]
        D_pseudo[end,end] -= inverse_weights[end]

    # Volume Integral for one Element, loop over all nodes
        for j in eachnode(dg), i in eachnode(dg)
            u_node = get_node_vars(u, equations, dg, i, j, element)
            q1_node = get_node_vars(q1, equations, dg, i, j, element)
            q2_node = get_node_vars(q2, equations, dg, i, j, element)

            flux1 = flux(u_node, 1, equations) - viscous_flux(u_node, q1_node, q2_node, 1, equations)
            flux2 = flux(u_node, 2, equations) - viscous_flux(u_node, q1_node, q2_node, 2, equations)

    # x direction
            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)

    # two point flux of the same node outside the loop. And Pseudo Strong Part for boundary nodes
            if i == 1 || i == nnodes(dg)
                # Here we can also use the phys. flux.
                contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
                integral_contribution = (2 * D_pseudo[i, i] + S[i, i]) *
                                        contravariant_flux1
                add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
            else
                contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
                integral_contribution = 2 * D_pseudo[i, i] * contravariant_flux1
                add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
            end

    # two point flux symmetry: estimate the volume flux only one time
            for ii in (i + 1):nnodes(dg)
                # arihtmetic mean for two point volume flux metrics
                Ja11_, Ja12_ = get_contravariant_vector(1, contravariant_vectors, ii, j, element)
                u_node_ii = get_node_vars(u, equations, dg, ii, j, element)
                q1_node_ii = get_node_vars(q1, equations, dg, ii, j, element)
                q2_node_ii = get_node_vars(q2, equations, dg, ii, j, element)

                flux1 = euler_and_viscous_flux(u_node, q1_node, q2_node, u_node_ii, q1_node_ii, q2_node_ii,
                                                1, equations, volume_integral.volume_flux)
                flux2 = euler_and_viscous_flux(u_node, q1_node, q2_node, u_node_ii, q1_node_ii, q2_node_ii,
                                                2, equations, volume_integral.volume_flux)
                contravariant_flux1 = 0.5 * (Ja11 + Ja11_) * flux1 + 0.5 * (Ja12 + Ja12_) * flux2
                integral_contribution = 2 * D_pseudo[i, ii] * contravariant_flux1
                add_to_node_vars!(du, integral_contribution, equations, dg, i,  j, element)
                integral_contribution = 2 * D_pseudo[ii, i] * contravariant_flux1
                add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
            end
    # y direction
            flux1 = flux(u_node, 1, equations) - viscous_flux(u_node, q1_node, q2_node, 1, equations)
            flux2 = flux(u_node, 2, equations) - viscous_flux(u_node, q1_node, q2_node, 2, equations)

            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)

    # two point flux of the same node outside the loop. And Pseudo Strong Part for boundary nodes
            if j == 1 || j == nnodes(dg)
                # Here we can also use the phys. flux.
                contra_flux = Ja22 * flux2 + Ja21 * flux1
                integral_contribution = (2 * D_pseudo[j, j] + S[j, j]) *
                                         contra_flux
                add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
            else
                contravariant_flux2 = Ja22 * flux2 + Ja21 * flux1
                integral_contribution = 2 * D_pseudo[j, j] * contravariant_flux2
                add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)
            end
    # two point flux symmetry: estimate the volume flux only one time
            for jj in (j + 1):nnodes(dg)
                # arithmetic mean for two point volume flux metrics
                Ja21_, Ja22_ = get_contravariant_vector(2, contravariant_vectors, i, jj, element)
                u_node_jj = get_node_vars(u, equations, dg, i, jj, element)
                q1_node_jj = get_node_vars(q1, equations, dg, i, jj, element)
                q2_node_jj = get_node_vars(q2, equations, dg, i, jj, element)

                flux1 = euler_and_viscous_flux(u_node, q1_node, q2_node, u_node_jj, q1_node_jj, q2_node_jj,
                                                    1, equations, volume_integral.volume_flux)
                flux2 = euler_and_viscous_flux(u_node, q1_node, q2_node, u_node_jj, q1_node_jj, q2_node_jj,
                                                    2, equations, volume_integral.volume_flux)
                contravariant_flux2 = 0.5 * (Ja22 + Ja22_) * flux2 + 0.5 * (Ja21 + Ja21_) * flux1
                integral_contribution =  2 * D_pseudo[j, jj] * contravariant_flux2
                add_to_node_vars!(du, integral_contribution, equations, dg, i, j,  element)
                integral_contribution =  2 * D_pseudo[jj, j] * contravariant_flux2
                add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
            end
        end
    end
end

function calc_interface_flux!(cache, u,
                              mesh::CurvedMesh{2},
                              equations, dg::DG)
    @unpack elements = cache

    @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2)"

    # Interfaces in x-direction (`orientation` = 1)
        calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[1, element],
                         element, 1, u, mesh, equations, dg, cache)

    # Interfaces in y-direction (`orientation` = 2)
        calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[2, element],
                         element, 2, u, mesh, equations, dg, cache)
    end

    return nothing
end

function calc_interface_flux_auxiliary!(nabla, cache, u,
    mesh::CurvedMesh{2},
    equations, dg::DG)
    @unpack elements = cache

    @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2)"

    # Interfaces in x-direction (`orientation` = 1)
        calc_interface_flux_auxiliary!(nabla, elements.surface_flux_values,
    elements.left_neighbors[1, element],
    element, 1, u, mesh, equations, dg, cache)

    # Interfaces in y-direction (`orientation` = 2)
        calc_interface_flux_auxiliary!(nabla, elements.surface_flux_values,
    elements.left_neighbors[2, element],
    element, 2, u, mesh, equations, dg, cache)
    end

    return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
                                      orientation, u,
                                      mesh::CurvedMesh{2}, equations,
                                      dg::DG, cache)
  # This is slow for LSA, but for some reason faster for Euler (see #519)
    if left_element <= 0 # left_element = 0 at boundaries
        return nothing
    end

    if typeof(equations) == CompressibleEulerEquations2D{Float64} && equations.viscous # BLZ: to add the central flux for viscous terms
        calc_interface_flux_viscous!(surface_flux_values, left_element, right_element,
                                      orientation, u,
                                      mesh, equations,
                                      dg, cache)
    else

        @unpack surface_flux = dg
        @unpack contravariant_vectors = cache.elements

        if typeof(equations) == AuxiliaryEquation # BLZ: to solve the auxiliary system
            surface_flux = flux_central_auxiliary
        end
        right_direction = 2 * orientation
        left_direction = right_direction - 1

        for i in eachnode(dg)
            if orientation == 1
                u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
                u_rr = get_node_vars(u, equations, dg, 1,          i, right_element)

      # First contravariant vector Ja^1 as SVector
                normal_vector = get_contravariant_vector(1, contravariant_vectors, 1, i, right_element)
            else # orientation == 2
                u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
                u_rr = get_node_vars(u, equations, dg, i, 1,          right_element)

      # Second contravariant vector Ja^2 as SVector
                normal_vector = get_contravariant_vector(2, contravariant_vectors, i, 1, right_element)
            end

            flux = surface_flux(u_ll, u_rr, normal_vector, equations)

            for v in eachvariable(equations)
                surface_flux_values[v, i, right_direction, left_element] = flux[v]
                surface_flux_values[v, i, left_direction, right_element] = flux[v]
            end
        end
    end

    return nothing
end

@inline function calc_interface_flux_viscous!(surface_flux_values, left_element, right_element,
    orientation, u,
    mesh::CurvedMesh{2}, equations,
    dg::DG, cache)
    # This is slow for LSA, but for some reason faster for Euler (see #519)
    if left_element <= 0 # left_element = 0 at boundaries
        return nothing
    end

    @unpack surface_flux = dg
    @unpack contravariant_vectors = cache.elements
    @unpack q1 = cache
    @unpack q2 = cache

    right_direction = 2 * orientation
    left_direction = right_direction - 1

    for i in eachnode(dg)
        if orientation == 1
            u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
            u_rr = get_node_vars(u, equations, dg, 1,          i, right_element)
            q1_ll = get_node_vars(q1, equations, dg, nnodes(dg), i, left_element)
            q1_rr = get_node_vars(q1, equations, dg, 1,          i, right_element)
            q2_ll = get_node_vars(q2, equations, dg, nnodes(dg), i, left_element)
            q2_rr = get_node_vars(q2, equations, dg, 1,          i, right_element)

    # First contravariant vector Ja^1 as SVector
            normal_vector = get_contravariant_vector(1, contravariant_vectors, 1, i, right_element)
        else # orientation == 2
            u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
            u_rr = get_node_vars(u, equations, dg, i, 1,          right_element)
            q1_ll = get_node_vars(q1, equations, dg, i, nnodes(dg), left_element)
            q1_rr = get_node_vars(q1, equations, dg, i, 1,          right_element)
            q2_ll = get_node_vars(q2, equations, dg, i, nnodes(dg), left_element)
            q2_rr = get_node_vars(q2, equations, dg, i, 1,          right_element)

    # Second contravariant vector Ja^2 as SVector
            normal_vector = get_contravariant_vector(2, contravariant_vectors, i, 1, right_element)
        end

        flux = surface_flux_plus_viscous_central(u_ll, u_rr, q1_ll, q1_rr, q2_ll, q2_rr, normal_vector, equations, surface_flux)

        for v in eachvariable(equations)
            surface_flux_values[v, i, right_direction, left_element] = flux[v]
            surface_flux_values[v, i, left_direction, right_element] = flux[v]
        end
    end

    return nothing
end

@inline function calc_interface_flux_auxiliary!(nabla, surface_flux_values, left_element, right_element,
    orientation, u,
    mesh::CurvedMesh{2}, equations,
    dg::DG, cache)
    # This is slow for LSA, but for some reason faster for Euler (see #519)
    if left_element <= 0 # left_element = 0 at boundaries
        return nothing
    end

    @unpack contravariant_vectors = cache.elements
    surface_flux = flux_central_auxiliary

    right_direction = 2 * orientation
    left_direction = right_direction - 1

    for i in eachnode(dg)
        if orientation == 1
            u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
            u_rr = get_node_vars(u, equations, dg, 1,          i, right_element)

    # First contravariant vector Ja^1 as SVector
            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, 1, i, right_element)
    # derivative in only one direction; y_η for u_x, -x_η for u_y, the other one is 0
            if nabla == 1 normal_vector = SVector(Ja11, 0) else normal_vector = SVector(0, Ja12) end
        else # orientation == 2
            u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
            u_rr = get_node_vars(u, equations, dg, i, 1,          right_element)

    # Second contravariant vector Ja^2 as SVector
            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, 1, right_element)
    # derivative in only one direction; -y_ξ for u_x, x_ξ for u_y, the other one is 0
            if nabla == 2 normal_vector = SVector(Ja22, 0) else normal_vector = normal_vector = SVector(0, Ja21) end
        end

        flux = surface_flux(u_ll, u_rr, normal_vector, equations)

        for v in eachvariable(equations)
            surface_flux_values[v, i, right_direction, left_element] = flux[v]
            surface_flux_values[v, i, left_direction, right_element] = flux[v]
        end
    end

    return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::CurvedMesh{2}, equations, dg::DG)
    @assert isperiodic(mesh)
end
function calc_boundary_flux_auxiliary!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
    mesh::CurvedMesh{2}, equations, dg::DG, nabla)
    @assert isperiodic(mesh)
    calc_boundary_flux_auxiliary!(cache, u, t, (boundary_condition, boundary_condition,
                        boundary_condition, boundary_condition), mesh, equations, dg, nabla)
end

function calc_boundary_flux!(cache, u, t, boundary_condition,
                             mesh::CurvedMesh{2}, equations, dg::DG)
    calc_boundary_flux!(cache, u, t,
                      (boundary_condition, boundary_condition,
                       boundary_condition, boundary_condition),
                      mesh, equations, dg)
end
function calc_boundary_flux_auxiliary!(cache, u, t, boundary_condition,
    mesh::CurvedMesh{2}, equations, dg::DG, nabla)
    calc_boundary_flux_auxiliary!(cache, u, t,
    (boundary_condition, boundary_condition,
    boundary_condition, boundary_condition),
    mesh, equations, dg, nabla)
end

function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple},
                             mesh::CurvedMesh{2}, equations, dg::DG)
    @unpack surface_flux = dg
    @unpack surface_flux_values = cache.elements
    linear_indices = LinearIndices(size(mesh))

    for cell_y in axes(mesh, 2)
    # Negative x-direction
        direction = 1
        element = linear_indices[begin, cell_y]

        for j in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction],
                                       mesh, equations, dg, cache,
                                       direction, (1, j), (j,), element)
        end

    # Positive x-direction
        direction = 2
        element = linear_indices[end, cell_y]

        for j in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction],
                                       mesh, equations, dg, cache,
                                       direction, (nnodes(dg), j), (j,), element)
        end
    end

    for cell_x in axes(mesh, 1)
    # Negative y-direction
        direction = 3
        element = linear_indices[cell_x, begin]

        for i in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction],
                                       mesh, equations, dg, cache,
                                       direction, (i, 1), (i,), element)
        end

    # Positive y-direction
        direction = 4
        element = linear_indices[cell_x, end]

        for i in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction],
                                       mesh, equations, dg, cache,
                                       direction, (i, nnodes(dg)), (i,), element)
        end
    end
end

# BLZ
function calc_boundary_flux_auxiliary!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple},
    mesh::CurvedMesh{2}, equations::AuxiliaryEquation, dg::DG, nabla)

    @unpack surface_flux_values = cache.elements
    linear_indices = LinearIndices(size(mesh))

    for cell_y in axes(mesh, 2)
    # Negative x-direction
        direction = 1
        element = linear_indices[begin, cell_y]
        # periodic boundaries computed as interface fluxes for the auxiliary equation
        if boundary_conditions[direction] == boundary_condition_periodic
            orientation = 1 # x axis
            left_element = element # in negative x direction the current element is the left element
            right_element = linear_indices[end, cell_y]
            calc_interface_flux_auxiliary!(nabla, surface_flux_values, left_element, right_element,
                                            orientation, u,
                                            mesh, equations,
                                            dg, cache)
        end

        for j in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
              boundary_conditions[direction],
              mesh, equations, dg, cache,
              direction, (1, j), (j,), element, nabla)
        end

    # Positive x-direction
        direction = 2
        element = linear_indices[end, cell_y]
        # periodic boundaries computed as interface fluxes for the auxiliary equation
        if boundary_conditions[direction] == boundary_condition_periodic
            orientation = 1 # x axis
            right_element = element # in positive x direction the current element is the right element
            left_element = linear_indices[begin, cell_y]
            calc_interface_flux_auxiliary!(nabla, surface_flux_values, left_element, right_element,
                                            orientation, u,
                                            mesh, equations,
                                            dg, cache)
        end

        for j in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
              boundary_conditions[direction],
              mesh, equations, dg, cache,
              direction, (nnodes(dg), j), (j,), element, nabla)
        end
    end

    for cell_x in axes(mesh, 1)
    # Negative y-direction
        direction = 3
        element = linear_indices[cell_x, begin]
        # periodic boundaries computed as interface fluxes for the auxiliary equation
        if boundary_conditions[direction] == boundary_condition_periodic
            orientation = 2 # y axis
            right_element = element # in negative y direction the current element is the right element
            left_element = linear_indices[cell_x, end]
            calc_interface_flux_auxiliary!(nabla, surface_flux_values, left_element, right_element,
                                            orientation, u,
                                            mesh, equations,
                                            dg, cache)
            # u_inner = get_node_vars(u, equations, dg, node_indices..., element)
            # flux = surface_flux(u_inner, uboundary, normal, equations)
            # for v in eachvariable(equations)
            #     surface_flux_values[v, i, right_direction, left_element] = flux[v]
            #     surface_flux_values[v, i, left_direction, right_element] = flux[v]
            # end
        end

        for i in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
              boundary_conditions[direction],
              mesh, equations, dg, cache,
              direction, (i, 1), (i,), element, nabla)
        end

    # Positive y-direction
        direction = 4
        element = linear_indices[cell_x, end]
        # periodic boundaries computed as interface fluxes for the auxiliary equation
        if boundary_conditions[direction] == boundary_condition_periodic
            orientation = 2 # y axis
            left_element = element # in positive y direction the current element is the left element
            right_element = linear_indices[cell_x, begin]
            calc_interface_flux_auxiliary!(nabla, surface_flux_values, left_element, right_element,
                                            orientation, u,
                                            mesh, equations,
                                            dg, cache)
        end

        for i in eachnode(dg)
            calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
              boundary_conditions[direction],
              mesh, equations, dg, cache,
              direction, (i, nnodes(dg)), (i,), element, nabla)
        end
    end
end

# BLZ modification neccessary for viscous equation, original function is in dg.jl
@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
    boundary_condition,
    mesh::CurvedMesh, equations, dg::DG, cache,
    direction, node_indices, surface_node_indices, element)
    @unpack node_coordinates, contravariant_vectors = cache.elements
    @unpack surface_flux = dg

    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
    x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)

  # Contravariant vector Ja^i is the normal vector
    normal = get_contravariant_vector(orientation, contravariant_vectors, node_indices..., element)
    if equations.viscous # boundaries including gradient needed for viscous equation
    @unpack q1 = cache
    @unpack q2 = cache
    q1_inner = get_node_vars(q1, equations, dg, node_indices..., element)
    q2_inner = get_node_vars(q2, equations, dg, node_indices..., element)
    flux = boundary_condition(u_inner, q1_inner, q2_inner, normal, direction, x, t, surface_flux, equations)
    else
        flux = boundary_condition(u_inner, normal, direction, x, t, surface_flux, equations)
    end
    for v in eachvariable(equations)
        surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
    end
end

@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
    boundary_condition,
    mesh::CurvedMesh, equations::AuxiliaryEquation, dg::DG, cache,
    direction, node_indices, surface_node_indices, element, nabla)
    @unpack node_coordinates, contravariant_vectors = cache.elements
    # central flux for gradient equation
    surface_flux = flux_central

    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
    x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)

  # Contravariant vector Ja^i is the normal vector
    normal = get_contravariant_vector(orientation, contravariant_vectors, node_indices..., element)

  # BLZ: derivative in only one direction, depending on orientation;
    # orientation = 2 : -y_ξ for u_x, x_ξ for u_y, the other one is 0
    # orientation = 1 : y_η for u_x, -x_η for u_y, the other one is 0
    if orientation == 2
        if nabla == 2
            normal_vector = SVector(normal[2], 0)
        else normal_vector = SVector(0, normal[1])
        end
    else # orientation == 1
        if nabla == 1
            normal_vector = SVector(normal[1], 0)
        else normal_vector = SVector(0, normal[2])
        end
    end
    # just set q1, q2 = 0 as placeholder
    flux = boundary_condition(u_inner, 0, 0, normal_vector, direction, x, t, surface_flux, equations)

    for v in eachvariable(equations)
        surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
    end
end

function apply_jacobian!(du,
                         mesh::Union{CurvedMesh{2},UnstructuredQuadMesh},
                         equations, dg::DG, cache)
    @unpack inverse_jacobian = cache.elements

    @threaded for element in eachelement(dg, cache)
        for j in eachnode(dg), i in eachnode(dg)
            factor = -inverse_jacobian[i, j, element]
            # factor = inverse_jacobian[i, j, element]

            for v in eachvariable(equations)
                du[v, i, j, element] *= factor
            end
        end
    end

    return nothing
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test

using IntervalSets

@testset "util" begin
    @testset "_car_cdr" begin
        @test BAT._car_cdr_impl() == ()
        @test BAT._car_cdr_impl(4, 3.5, -1) == (4, (3.5, -1))
        @test BAT._car_cdr((3, -8, 2)) == (3, (-8, 2))
    end

    @testset "_all_lteq" begin
        @test BAT._all_lteq(-0.4, -0.4:0.4:0.8 ,0.8)
        @test !BAT._all_lteq(-0.4, -0.4:0.4:0.8 ,0.7)
        @test !BAT._all_lteq(-0.3, -0.4:0.4:0.8 ,0.8)

        @test BAT._all_lteq([0.0, 0.5, -0.4], [0.3, 0.5, -0.4] , [0.3, 0.6, -0.3])
        @test !BAT._all_lteq([0.5, -0.4], [0.49, -0.4] , [0.6, -0.3])
        @test_throws DimensionMismatch BAT._all_lteq(
            [0.5, -0.4], [0.3, 0.5, -0.4] , [0.3, 0.6, -0.3])

        @test BAT._all_in_ui(ones(2, 2))
        @test !BAT._all_in_ui([0.4, 0.2, 1.1])
    end

    @testset "fromui" begin
        @test fromui(0.5, -1.0, 1.0) ≈ 0.0
        @test fromui(0.5, -0.5, 5.0) ≈ 2.25
        @test fromui(0.1, -0.5, 1.0) ≈ -0.35
        @test_throws ArgumentError fromui(-0.1, -1.0, 1.0)
        @test_throws ArgumentError fromui(1.1, -0.5, 1.0)

        @test fromui(0.5, ClosedInterval(-0.5, 1.0)) ≈ 0.25

        inv_fromui = @inferred inv(fromui)
        @test inv_fromui  == BAT.inv_fromui
        @test inv(inv_fromui) == fromui

        @test inv_fromui(2.25, -0.5, 5.0) ≈ 0.5

        @test inv_fromui(-0.35, -0.5, 1.0) ≈ 0.1
        @test inv_fromui(5.0, -0.5, 5.0) ≈ 1.0
        @test inv_fromui(-0.5, -0.5, 5.0) ≈ 0.0
        @test_throws ArgumentError inv_fromui(-0.2, -0.1, 3)
        @test_throws ArgumentError inv_fromui(3.1, -0.1, 3)

        @test inv_fromui(-0.35, ClosedInterval(-0.5, 1.0)) ≈ 0.1
    end
end


text:
function crosscorr_input(df::DataFrame, ycol::Symbol, features::Array{Symbol, 1}, lags)
    output_dir = "/mnt/crosscorr/"
    Base.Filesystem.mkpath(output_dir)

    for fea in features
        cor = crosscor(df[:, ycol], df[:, fea], lags)
        plot_corr_OU(lags, cor, "$(string(ycol)) & $(string(fea))", output_dir, "crosscor_$(string(ycol))_$(string(fea))")
    end
end

function autocorr_input(df::DataFrame, ycol::Symbol, lags)
    output_dir = "/mnt/autocorr/"
    Base.Filesystem.mkpath(output_dir)

    cor = autocor(df[:, ycol], lags)
    plot_corr_OU(cor, "$(string(ycol))", output_dir, "autocor_$(string(ycol))")

end


text:
struct PeriodicBoundaryCondition <: AbstractBoundaryCondition end
const Periodic = PeriodicBoundaryCondition()