text:
abstract type AbstractSpectralStyle <: Fields.AbstractFieldStyle end

struct SpectralStyle <: AbstractSpectralStyle end

struct CompositeSpectralStyle <: AbstractSpectralStyle end

"""
    SpectralElementOperator

Represents an operation that is applied to each element.

Subtypes `Op` of this should define the following:
- [`operator_return_eltype(::Op, ElTypes...)`](@ref)
- [`allocate_work(::Op, args...)`](@ref)
- [`apply_slab(::Op, work, args...)`](@ref)

Additionally, the result type `OpResult <: OperatorSlabResult` of `apply_slab` should define `get_node(::OpResult, i, j)`.
"""
abstract type SpectralElementOperator end

"""
    operator_axes(space)

Return a tuple of the axis indicies a given field operator works over.
"""
function operator_axes end

operator_axes(space::Spaces.AbstractSpace) = ()
operator_axes(space::Spaces.SpectralElementSpace1D) = (1,)
operator_axes(space::Spaces.SpectralElementSpace2D) = (1, 2)
operator_axes(space::Spaces.ExtrudedFiniteDifferenceSpace) =
    operator_axes(space.horizontal_space)

"""
    SpectralBroadcasted{Style}(op, args[,axes[, work]])

This is similar to a `Base.Broadcast.Broadcasted` object, except it contains space for an intermediate `work` storage.

This is returned by `Base.Broadcast.broadcasted(op::SpectralElementOperator)`.
"""
struct SpectralBroadcasted{Style, Op, Args, Axes, InputSpace} <:
       Base.AbstractBroadcasted
    op::Op
    args::Args
    axes::Axes
    input_space::InputSpace
end
SpectralBroadcasted{Style}(
    op::Op,
    args::Args,
    axes::Axes = nothing,
    input_space::InputSpace = nothing,
) where {Style, Op, Args, Axes, InputSpace} =
    SpectralBroadcasted{Style, Op, Args, Axes, InputSpace}(
        op,
        args,
        axes,
        input_space,
    )

input_space(arg) = axes(arg)
input_space(::SpectralElementOperator, space) = space

input_space(sbc::SpectralBroadcasted) =
    isnothing(sbc.input_space) ? input_space(sbc.op, map(axes, sbc.args)...) :
    sbc.input_space

return_space(::SpectralElementOperator, space) = space

Base.axes(sbc::SpectralBroadcasted) =
    isnothing(sbc.axes) ? return_space(sbc.op, map(axes, sbc.args)...) :
    sbc.axes

Base.Broadcast.broadcasted(op::SpectralElementOperator, args...) =
    Base.Broadcast.broadcasted(SpectralStyle(), op, args...)

Base.Broadcast.broadcasted(
    ::SpectralStyle,
    op::SpectralElementOperator,
    args...,
) = SpectralBroadcasted{SpectralStyle}(op, args)

Base.eltype(sbc::SpectralBroadcasted) =
    operator_return_eltype(sbc.op, map(eltype, sbc.args)...)

function Base.Broadcast.instantiate(
    sbc::SpectralBroadcasted{Style},
) where {Style}
    op = sbc.op
    # recursively instantiate the arguments to allocate intermediate work arrays
    args = map(Base.Broadcast.instantiate, sbc.args)
    # axes: same logic as Broadcasted
    if sbc.axes isa Nothing # Not done via dispatch to make it easier to extend instantiate(::Broadcasted{Style})
        axes = Base.axes(sbc)
    else
        axes = sbc.axes
        Base.Broadcast.check_broadcast_axes(axes, args...)
    end
    if sbc.input_space isa Nothing
        inspace = input_space(sbc)
    else
        inspace = sbc.input_space
    end
    op = typeof(op)(axes)
    return SpectralBroadcasted{Style}(op, args, axes, inspace)
end

function Base.Broadcast.instantiate(
    bc::Base.Broadcast.Broadcasted{Style},
) where {Style <: AbstractSpectralStyle}
    # recursively instantiate the arguments to allocate intermediate work arrays
    args = map(Base.Broadcast.instantiate, bc.args)
    # axes: same logic as Broadcasted
    if bc.axes isa Nothing # Not done via dispatch to make it easier to extend instantiate(::Broadcasted{Style})
        axes = Base.Broadcast.combine_axes(args...)
    else
        axes = bc.axes
        Base.Broadcast.check_broadcast_axes(axes, args...)
    end
    return Base.Broadcast.Broadcasted{Style}(bc.f, args, axes)
end

function Base.similar(sbc::SpectralBroadcasted, ::Type{Eltype}) where {Eltype}
    space = axes(sbc)
    return Field(Eltype, space)
end

function Base.similar(
    bc::Base.Broadcast.Broadcasted{<:AbstractSpectralStyle},
    ::Type{Eltype},
) where {Eltype}
    space = axes(bc)
    return Field(Eltype, space)
end

function Base.copy(sbc::SpectralBroadcasted)
    # figure out return type
    dest = similar(sbc, eltype(sbc))
    # allocate dest
    copyto!(dest, sbc)
end
Base.Broadcast.broadcastable(sbc::SpectralBroadcasted) = sbc

function Base.Broadcast.materialize(sbc::SpectralBroadcasted)
    copy(Base.Broadcast.instantiate(sbc))
end

function Base.copyto!(field_out::Field, sbc::SpectralBroadcasted)
    data_out = Fields.field_values(field_out)
    space = axes(field_out)
    Nv = Spaces.nlevels(space)
    Nh = Topologies.nlocalelems(Spaces.topology(space))
    for h in 1:Nh, v in 1:Nv
        slab_out = slab(field_out, v, h)
        out_slab_space = slab(axes(sbc), v, h)
        in_slab_space = slab(input_space(sbc), v, h)
        slab_args = map(arg -> _apply_slab(slab(arg, v, h), v, h), sbc.args)
        copy_slab!(
            slab_out,
            apply_slab(sbc.op, out_slab_space, in_slab_space, slab_args...),
        )
    end
    return field_out
end

function Base.Broadcast.materialize!(dest, sbc::SpectralBroadcasted)
    copyto!(dest, Base.Broadcast.instantiate(sbc))
end

function slab(
    sbc::SpectralBroadcasted{Style},
    inds...,
) where {Style <: SpectralStyle}
    _args = map(a -> slab(a, inds...), sbc.args)
    _axes = slab(axes(sbc), inds...)
    SpectralBroadcasted{Style}(sbc.op, _args, _axes, sbc.input_space)
end

function slab(
    bc::Base.Broadcast.Broadcasted{Style},
    inds...,
) where {Style <: AbstractSpectralStyle}
    _args = map(a -> slab(a, inds...), bc.args)
    _axes = slab(axes(bc), inds...)
    Base.Broadcast.Broadcasted{Style}(bc.f, _args, _axes)
end

function Base.copyto!(
    field_out::Field,
    bc::Base.Broadcast.Broadcasted{Style},
) where {Style <: AbstractSpectralStyle}
    data_out = Fields.field_values(field_out)
    space = axes(field_out)
    Nv = Spaces.nlevels(space)
    Nh = Topologies.nlocalelems(Spaces.topology(space))
    for h in 1:Nh, v in 1:Nv
        slab_out = slab(field_out, v, h)
        slab_args = map(arg -> _apply_slab(slab(arg, v, h), v, h), bc.args)
        copy_slab!(
            slab_out,
            Base.Broadcast.Broadcasted{Style}(bc.f, slab_args, axes(slab_out)),
        )
    end
    return field_out
end

function copy_slab!(slab_out::Fields.SlabField1D, res)
    space = axes(slab_out)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)
    @inbounds for i in 1:Nq
        set_node!(slab_out, i, get_node(res, i))
    end
    return slab_out
end

function copy_slab!(slab_out::Fields.SlabField2D, res)
    space = axes(slab_out)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)
    @inbounds for j in 1:Nq, i in 1:Nq
        set_node!(slab_out, i, j, get_node(res, i, j))
    end
    return slab_out
end

# 1D get/set node

# 1D intermediate slab data types
Base.Base.@propagate_inbounds function get_node(v::MVector, i)
    v[i]
end

Base.Base.@propagate_inbounds function get_node(v::SVector, i)
    v[i]
end

Base.Base.@propagate_inbounds function get_node(scalar, i)
    scalar[]
end

Base.@propagate_inbounds function get_node(field::Fields.SlabField1D, i)
    getindex(Fields.field_values(field), i)
end

@inline function get_node(bc::Base.Broadcast.Broadcasted, i)
    args = map(arg -> get_node(arg, i), bc.args)
    bc.f(args...)
end

Base.@propagate_inbounds set_node!(field::Fields.SlabField1D, i, val) =
    setindex!(Fields.field_values(field), val, i)

# 2D get/set node
Base.@propagate_inbounds function get_node(v::MMatrix, i, j)
    v[i, j]
end

Base.@propagate_inbounds function get_node(v::SMatrix, i, j)
    v[i, j]
end

Base.@propagate_inbounds function get_node(scalar, i, j)
    scalar[]
end

Base.@propagate_inbounds function get_node(field::Fields.SlabField2D, i, j)
    getindex(Fields.field_values(field), i, j)
end

@inline function get_node(bc::Base.Broadcast.Broadcasted, i, j)
    args = map(arg -> get_node(arg, i, j), bc.args)
    bc.f(args...)
end

Base.@propagate_inbounds function set_node!(
    field::Fields.SlabField2D,
    i,
    j,
    val,
)
    setindex!(Fields.field_values(field), val, i, j)
end

# Broadcast recursive machinery
_apply_slab(x, inds...) = x

_apply_slab(sbc::SpectralBroadcasted, inds...) = apply_slab(
    sbc.op,
    slab(axes(sbc), inds...),
    slab(input_space(sbc), inds...),
    map(a -> _apply_slab(a, inds...), sbc.args)...,
)

_apply_slab(bc::Base.Broadcast.Broadcasted{CompositeSpectralStyle}, inds...) =
    Base.Broadcast.Broadcasted{CompositeSpectralStyle}(
        bc.f,
        map(a -> _apply_slab(a, inds...), bc.args),
        bc.axes,
    )

function Base.Broadcast.BroadcastStyle(
    ::Type{SB},
) where {SB <: SpectralBroadcasted}
    CompositeSpectralStyle()
end

function Base.Broadcast.BroadcastStyle(
    ::CompositeSpectralStyle,
    ::Fields.AbstractFieldStyle,
)
    CompositeSpectralStyle()
end

"""
    Divergence()

Computes the "strong" divergence of a vector field `v`.

We compute the divergence as

    [∂(Jv¹)/∂ξ¹ + ∂(Jv²)/∂ξ²] / J

where `J` is the Jacobian determinant, `vⁱ` is the `i`th contravariant component of `v`.

This is discretized at the quadrature points as

    I{[∂(I{Jv¹})/∂ξ¹ + ∂(I{Jv²})/∂ξ²] / J}

where `I{x}` is the interpolation operator applied to a field `x`.

## References
 - Taylor and Fournier (2010), equation 15
"""
struct Divergence{I} <: SpectralElementOperator end
Divergence() = Divergence{()}()
Divergence{()}(space) = Divergence{operator_axes(space)}()

operator_return_eltype(op::Divergence, S) =
    RecursiveApply.rmaptype(Geometry.divergence_result_type, S)

function apply_slab(op::Divergence{(1,)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        Jv¹ =
            local_geometry.J ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant1(v, local_geometry),
                get_node(slab_data, i),
            )
        for ii in 1:Nq
            out[ii] = out[ii] ⊞ (D[ii, i] ⊠ Jv¹)
        end
    end
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        out[i] = RecursiveApply.rdiv(out[i], local_geometry.J)
    end
    return SVector(out)
end

function apply_slab(op::Divergence{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        Jv¹ =
            local_geometry.J ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant1(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for ii in 1:Nq
            out[ii, j] = out[ii, j] ⊞ (D[ii, i] ⊠ Jv¹)
        end
        Jv² =
            local_geometry.J ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant2(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for jj in 1:Nq
            out[i, jj] = out[i, jj] ⊞ (D[jj, j] ⊠ Jv²)
        end
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.J)
    end
    return SMatrix(out)
end

"""
    WeakDivergence()

Computes the divergence of a vector field `v` weakly.

This spolves the variational problem of finding `θ` such that

    ⟨ϕ, J θ⟩ = - ⟨∂ϕ/∂ξ¹, J v¹⟩ + ⟨∂ϕ/∂ξ², J v²⟩

for all `ϕ` (which arises by integration by parts).

Discretely it is equivalent to

    - J \\ (D₁' * W * J * v¹ + D₂' * W * J * v²)

where

 - `J` is the diagonal Jacobian matrix
 - `W` is the diagonal matrix of quadrature weights
 - `D₁` and `D₂` are the discrete derivative matrices along the first and second dimensions.
"""
struct WeakDivergence{I} <: SpectralElementOperator end
WeakDivergence() = WeakDivergence{()}()
WeakDivergence{()}(space) = WeakDivergence{operator_axes(space)}()

operator_return_eltype(op::WeakDivergence, S) =
    RecursiveApply.rmaptype(Geometry.divergence_result_type, S)

function apply_slab(op::WeakDivergence{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    RT = operator_return_eltype(op, eltype(slab_data))
    # allocate temp output
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        WJv¹ =
            local_geometry.WJ ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant1(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for ii in 1:Nq
            out[ii, j] = out[ii, j] ⊞ (D[i, ii] ⊠ WJv¹)
        end
        WJv² =
            local_geometry.WJ ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant2(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for jj in 1:Nq
            out[i, jj] = out[i, jj] ⊞ (D[j, jj] ⊠ WJv²)
        end
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], ⊟(local_geometry.WJ))
    end
    return SMatrix(out)
end

"""
    Gradient()

Compute the (strong) gradient on each element via the chain rule:

    ∂f/∂xⁱ = ∂f/∂ξʲ * ∂ξʲ/∂xⁱ
"""
struct Gradient{I} <: SpectralElementOperator end
Gradient() = Gradient{()}()
Gradient{()}(space) = Gradient{operator_axes(space)}()

operator_return_eltype(::Gradient{I}, S) where {I} =
    RecursiveApply.rmaptype(T -> Geometry.gradient_result_type(Val(I), T), S)

function apply_slab(op::Gradient{(1,)}, slab_space, _, slab_data)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for i in 1:Nq
        x = get_node(slab_data, i)
        for ii in 1:Nq
            ∂f∂ξ = Geometry.Covariant1Vector(D[ii, i]) ⊗ x
            out[ii] += ∂f∂ξ
        end
    end
    return SVector(out)
end

function apply_slab(op::Gradient{(1, 2)}, slab_space, _, slab_data)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        x = get_node(slab_data, i, j)
        for ii in 1:Nq
            ∂f∂ξ₁ = Geometry.Covariant12Vector(D[ii, i], zero(eltype(D))) ⊗ x
            out[ii, j] = out[ii, j] ⊞ ∂f∂ξ₁
        end
        for jj in 1:Nq
            ∂f∂ξ₂ = Geometry.Covariant12Vector(zero(eltype(D)), D[jj, j]) ⊗ x
            out[i, jj] = out[i, jj] ⊞ ∂f∂ξ₂
        end
    end
    return SMatrix(out)
end

"""
    WeakGradient()

Compute the (strong) gradient on each element via the chain rule:

    ∂f/∂xⁱ = ∂f/∂ξʲ * ∂ξʲ/∂xⁱ
"""
struct WeakGradient{I} <: SpectralElementOperator end
WeakGradient() = WeakGradient{()}()
WeakGradient{()}(space) = WeakGradient{operator_axes(space)}()

operator_return_eltype(::WeakGradient{I}, S) where {I} =
    RecursiveApply.rmaptype(T -> Geometry.gradient_result_type(Val(I), T), S)

function apply_slab(op::WeakGradient{(1,)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        W = local_geometry.WJ / local_geometry.J
        Wx = W ⊠ get_node(slab_data, i)
        for ii in 1:Nq
            Dᵀ₁Wf = Geometry.Covariant1Vector(D[i, ii]) ⊗ Wx
            out[ii] = out[ii] ⊟ Dᵀ₁Wf
        end
    end
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        W = local_geometry.WJ / local_geometry.J
        out[i] = RecursiveApply.rdiv(out[i], W)
    end
    return SVector(out)
end

function apply_slab(op::WeakGradient{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        W = local_geometry.WJ / local_geometry.J
        Wx = W ⊠ get_node(slab_data, i, j)
        for ii in 1:Nq
            Dᵀ₁Wf = Geometry.Covariant12Vector(D[i, ii], zero(eltype(D))) ⊗ Wx
            out[ii, j] = out[ii, j] ⊟ Dᵀ₁Wf
        end
        for jj in 1:Nq
            Dᵀ₂Wf = Geometry.Covariant12Vector(zero(eltype(D)), D[j, jj]) ⊗ Wx
            out[i, jj] = out[i, jj] ⊟ Dᵀ₂Wf
        end
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        W = local_geometry.WJ / local_geometry.J
        out[i, j] = RecursiveApply.rdiv(out[i, j], W)
    end
    return SMatrix(out)
end

abstract type CurlSpectralElementOperator <: SpectralElementOperator end

"""
    Curl()

Compute the (strong) curl on each element
"""
struct Curl{I} <: CurlSpectralElementOperator end
Curl() = Curl{()}()
Curl{()}(space) = Curl{operator_axes(space)}()

operator_return_eltype(::Curl{(1, 2)}, S) =
    RecursiveApply.rmaptype(T -> Geometry.curl_result_type(T), S)

function apply_slab(op::Curl{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    if RT <: Geometry.Contravariant3Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            v₁ = Geometry.covariant1(
                get_node(slab_data, i, j),
                slab_local_geometry[i, j],
            )
            for jj in 1:Nq
                D₂v₁ = D[jj, j] ⊠ v₁
                out[i, jj] =
                    out[i, jj] ⊞ Geometry.Contravariant3Vector(
                        ⊟(D₂v₁),
                        slab_local_geometry[i, jj],
                    )
            end
            v₂ = Geometry.covariant2(
                get_node(slab_data, i, j),
                slab_local_geometry[i, j],
            )
            for ii in 1:Nq
                D₁v₂ = D[ii, i] ⊠ v₂
                out[ii, j] =
                    out[ii, j] ⊞ Geometry.Contravariant3Vector(
                        D₁v₂,
                        slab_local_geometry[ii, j],
                    )
            end
        end
    elseif RT <: Geometry.Contravariant12Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            v₃ = Geometry.covariant3(
                get_node(slab_data, i, j),
                slab_local_geometry[i, j],
            )
            for ii in 1:Nq
                D₁v₃ = D[ii, i] ⊠ v₃
                out[ii, j] =
                    out[ii, j] ⊞
                    Geometry.Contravariant12Vector(zero(D₁v₃), ⊟(D₁v₃))
            end
            for jj in 1:Nq
                D₂v₃ = D[jj, j] ⊠ v₃
                out[i, jj] =
                    out[i, jj] ⊞
                    Geometry.Contravariant12Vector(D₂v₃, zero(D₂v₃))
            end
        end
    else
        error("invalid return type: $RT")
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.J)
    end
    return SMatrix(out)
end

"""
    WeakCurl()

Compute the weak curl on each element
"""
struct WeakCurl{I} <: CurlSpectralElementOperator end
WeakCurl() = WeakCurl{()}()
WeakCurl{()}(space) = WeakCurl{operator_axes(space)}()

operator_return_eltype(::WeakCurl{(1, 2)}, S) =
    RecursiveApply.rmaptype(T -> Geometry.curl_result_type(T), S)

function apply_slab(op::WeakCurl{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    if RT <: Geometry.Contravariant3Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            local_geometry = slab_local_geometry[i, j]
            W = local_geometry.WJ / local_geometry.J
            Wv₁ =
                W ⊠ Geometry.covariant1(
                    get_node(slab_data, i, j),
                    slab_local_geometry[i, j],
                )
            for jj in 1:Nq
                Dᵀ₂Wv₁ = D[j, jj] ⊠ Wv₁
                out[i, jj] =
                    out[i, jj] ⊞ Geometry.Contravariant3Vector(
                        Dᵀ₂Wv₁,
                        slab_local_geometry[i, jj],
                    )
            end
            Wv₂ =
                W ⊠ Geometry.covariant2(
                    get_node(slab_data, i, j),
                    slab_local_geometry[i, j],
                )
            for ii in 1:Nq
                Dᵀ₁Wv₂ = D[i, ii] ⊠ Wv₂
                out[ii, j] =
                    out[ii, j] ⊞ Geometry.Contravariant3Vector(
                        ⊟(Dᵀ₁Wv₂),
                        slab_local_geometry[ii, j],
                    )
            end
        end
    elseif RT <: Geometry.Contravariant12Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            local_geometry = slab_local_geometry[i, j]
            W = local_geometry.WJ / local_geometry.J
            Wv₃ =
                W ⊠ Geometry.covariant3(
                    get_node(slab_data, i, j),
                    slab_local_geometry[i, j],
                )
            for ii in 1:Nq
                Dᵀ₁Wv₃ = D[i, ii] ⊠ Wv₃
                out[ii, j] =
                    out[ii, j] ⊞
                    Geometry.Contravariant12Vector(zero(Dᵀ₁Wv₃), Dᵀ₁Wv₃)
            end
            for jj in 1:Nq
                Dᵀ₂Wv₃ = D[j, jj] ⊠ Wv₃
                out[i, jj] =
                    out[i, jj] ⊞
                    Geometry.Contravariant12Vector(⊟(Dᵀ₂Wv₃), zero(Dᵀ₂Wv₃))
            end
        end
    else
        error("invalid return type: $RT")
    end
    for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.WJ)
    end
    return SMatrix(out)
end

# interplation / restriction
abstract type TensorOperator <: SpectralElementOperator end

input_space(op::TensorOperator, inspace) = inspace
return_space(op::TensorOperator, inspace) = op.space
operator_return_eltype(op::TensorOperator, S) = S

"""
    Interpolate(space)

Computes the projection of a field to a higher degree polynomial space. `space`
must be on the same element mesh as the field, but have equal or higher
polynomial degree.

    ⟨ϕ, J θ⟩ = ⟨ϕ, J σ⟩

where `ϕ` and `θ` are on the higher degree space.

Discretely it is equivalent to

    I σ

where `I` is the interpolation matrix.
"""
struct Interpolate{I, S} <: TensorOperator
    space::S
end
Interpolate(space) = Interpolate{operator_axes(space), typeof(space)}(space)

function apply_slab(
    op::Interpolate{(1,)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    Imat = Quadratures.interpolation_matrix(FT, QS_out, QS_in)
    S = eltype(slab_data)
    slab_data_out = MVector{Nq_out, S}(undef)
    @inbounds for i in 1:Nq_out
        # manually inlined rmatmul with slab_getnode
        r = Imat[i, 1] ⊠ get_node(slab_data, 1)
        for ii in 2:Nq_in
            r = RecursiveApply.rmuladd(Imat[i, ii], get_node(slab_data, ii), r)
        end
        slab_data_out[i] = r
    end
    return slab_data_out
end

function apply_slab(
    op::Interpolate{(1, 2)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    Imat = Quadratures.interpolation_matrix(FT, QS_out, QS_in)
    S = eltype(slab_data)
    # temporary storage
    temp = MArray{Tuple{Nq_out, Nq_in}, S, 2, Nq_out * Nq_in}(undef)
    slab_data_out = MArray{Tuple{Nq_out, Nq_out}, S, 2, Nq_out * Nq_out}(undef)
    @inbounds for j in 1:Nq_in, i in 1:Nq_out
        # manually inlined rmatmul1 with slab get_node
        # we do this to remove one allocated intermediate array
        r = Imat[i, 1] ⊠ get_node(slab_data, 1, j)
        for ii in 2:Nq_in
            r = RecursiveApply.rmuladd(
                Imat[i, ii],
                get_node(slab_data, ii, j),
                r,
            )
        end
        temp[i, j] = r
    end
    @inbounds for j in 1:Nq_out, i in 1:Nq_out
        slab_data_out[i, j] = RecursiveApply.rmatmul2(Imat, temp, i, j)
    end
    return SMatrix(slab_data_out)
end

"""
    Restrict(space)

Computes the projection of a field to a lower degree polynomial space. `space`
must be on the same element mesh as the field, but have lower polynomial degree.

    ⟨ϕ, J θ⟩ = ⟨ϕ, J σ⟩

where `ϕ` and `θ` are on the lower degree space.

Discretely it is equivalent to

    (JWr) \\  I' (JW) σ

where `I` is the interpolation matrix, and `JWr` is the Jacobian multiplied by
quadrature weights on the lower-degree space.
"""
struct Restrict{I, S} <: TensorOperator
    space::S
end
Restrict(space) = Restrict{operator_axes(space), typeof(space)}(space)

function apply_slab(
    op::Restrict{(1,)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    ImatT = Quadratures.interpolation_matrix(FT, QS_in, QS_out)' # transpose
    S = eltype(slab_data)
    slab_data_out = MVector{Nq_out, S}(undef)
    slab_local_geometry_in = Spaces.local_geometry_data(slab_space_in)
    slab_local_geometry_out = Spaces.local_geometry_data(slab_space_out)
    WJ_in = slab_local_geometry_in.WJ
    WJ_out = slab_local_geometry_out.WJ
    @inbounds for i in 1:Nq_out
        # manually inlined rmatmul with slab get_node
        r = ImatT[i, 1] ⊠ (WJ_in[1] ⊠ get_node(slab_data, 1))
        for ii in 2:Nq_in
            WJ_node = WJ_in[ii] ⊠ get_node(slab_data, ii)
            r = RecursiveApply.rmuladd(ImatT[i, ii], WJ_node, r)
        end
        slab_data_out[i] = RecursiveApply.rdiv(r, WJ_out[i])
    end
    return slab_data_out
end

function apply_slab(
    op::Restrict{(1, 2)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    ImatT = Quadratures.interpolation_matrix(FT, QS_in, QS_out)' # transpose
    S = eltype(slab_data)
    # temporary storage
    temp = MArray{Tuple{Nq_out, Nq_in}, S, 2, Nq_out * Nq_in}(undef)
    slab_data_out = MArray{Tuple{Nq_out, Nq_out}, S, 2, Nq_out * Nq_out}(undef)
    slab_local_geometry_in = Spaces.local_geometry_data(slab_space_in)
    slab_local_geometry_out = Spaces.local_geometry_data(slab_space_out)
    WJ_in = slab_local_geometry_in.WJ
    WJ_out = slab_local_geometry_out.WJ
    @inbounds for j in 1:Nq_in, i in 1:Nq_out
        # manually inlined rmatmul1 with slab get_node
        r = ImatT[i, 1] ⊠ (WJ_in[1, j] ⊠ get_node(slab_data, 1, j))
        for ii in 2:Nq_in
            WJ_node = WJ_in[ii, j] ⊠ get_node(slab_data, ii, j)
            r = RecursiveApply.rmuladd(ImatT[i, ii], WJ_node, r)
        end
        temp[i, j] = r
    end
    @inbounds for j in 1:Nq_out, i in 1:Nq_out
        slab_data_out[i, j] = RecursiveApply.rdiv(
            RecursiveApply.rmatmul2(ImatT, temp, i, j),
            WJ_out[i, j],
        )
    end
    return SMatrix(slab_data_out)
end

"""
    tensor_product!(out, in, M)
    tensor_product!(inout, M)

Computes the tensor product `out = (M ⊗ M) * in` on each element.
"""
function tensor_product! end

function tensor_product!(
    out::DataLayouts.Data1DX{S, Ni_out},
    in::DataLayouts.Data1DX{S, Ni_in},
    M::SMatrix{Ni_out, Ni_in},
) where {S, Ni_out, Ni_in}
    (_, _, _, Nv_in, Nh_in) = size(in)
    (_, _, _, Nv_out, Nh_out) = size(out)
    # TODO: assumes the same number of levels (horizontal only)
    @assert Nv_in == Nv_out
    @assert Nh_in == Nh_out
    for h in 1:Nh_out, v in 1:Nv_out
        in_slab = slab(in, v, h)
        out_slab = slab(out, v, h)
        @inbounds for i in 1:Ni_out
            r = M[i, 1] ⊠ in_slab[1]
            for ii in 2:Ni_in
                r = RecursiveApply.rmuladd(M[i, ii], in_slab[ii], r)
            end
            out_slab[i] = r
        end
    end
    return out
end

function tensor_product!(
    out::DataLayouts.Data2D{S, Nij_out},
    in::DataLayouts.Data2D{S, Nij_in},
    M::SMatrix{Nij_out, Nij_in},
) where {S, Nij_out, Nij_in}

    Nh = length(in)
    @assert Nh == length(out)

    # temporary storage
    temp = MArray{Tuple{Nij_out, Nij_in}, S, 2, Nij_out * Nij_in}(undef)

    for h in 1:Nh
        in_slab = slab(in, h)
        out_slab = slab(out, h)
        for j in 1:Nij_in, i in 1:Nij_out
            temp[i, j] = RecursiveApply.rmatmul1(M, in_slab, i, j)
        end
        for j in 1:Nij_out, i in 1:Nij_out
            out_slab[i, j] = RecursiveApply.rmatmul2(M, temp, i, j)
        end
    end
    return out
end

"""
    tensor_product!(out, in, M)
    tensor_product!(inout, M)

Computes the tensor product `out = (M ⊗ M) * in` on each element.
"""
function tensor_product!(
    out_slab::DataLayouts.DataSlab2D{S, Nij_out},
    in_slab::DataLayouts.DataSlab2D{S, Nij_in},
    M::SMatrix{Nij_out, Nij_in},
) where {S, Nij_out, Nij_in}
    # temporary storage
    temp = MArray{Tuple{Nij_out, Nij_in}, S, 2, Nij_out * Nij_in}(undef)
    for j in 1:Nij_in, i in 1:Nij_out
        temp[i, j] = RecursiveApply.rmatmul1(M, in_slab, i, j)
    end
    for j in 1:Nij_out, i in 1:Nij_out
        out_slab[i, j] = RecursiveApply.rmatmul2(M, temp, i, j)
    end
    return out_slab
end

function tensor_product!(
    inout::Data2D{S, Nij},
    M::SMatrix{Nij, Nij},
) where {S, Nij}
    tensor_product!(inout, inout, M)
end

"""
    slab_gradient!(∇data, data, space)

Compute the gradient on each element via the chain rule:

    ∂f/∂xⁱ = ∂f/∂ξʲ * ∂ξʲ/∂xⁱ
"""
function slab_gradient!(∇data, data, space)
    # all derivatives calculated in the reference local geometry FT precision
    FT = Spaces.undertype(space)
    D = Quadratures.differentiation_matrix(FT, space.quadrature_style)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)

    # for each element in the element stack
    Nh = length(data)
    for h in 1:Nh
        ∇data_slab = slab(∇data, h)
        # TODO: can we just call materialize(slab(data,h))
        data_slab = slab(data, h)
        local_geometry_slab = slab(space.local_geometry, h)

        for i in 1:Nq, j in 1:Nq
            # TODO: materialize data_slab once
            # on GPU this would be done to shared memory, then synchronize threads
            local_geometry = local_geometry_slab[i, j]

            # compute covariant derivatives
            ∂f∂ξ₁ = RecursiveApply.rmatmul1(D, data_slab, i, j)
            ∂f∂ξ₂ = RecursiveApply.rmatmul2(D, data_slab, i, j)
            ∂f∂ξ = RecursiveApply.rmap(Covariant12Vector, ∂f∂ξ₁, ∂f∂ξ₂)

            # convert to desired basis
            ∇data_slab[i, j] = RecursiveApply.rmap(
                x -> Cartesian12Vector(x, local_geometry),
                ∂f∂ξ,
            )
        end
    end
    return ∇data
end

"""
    slab_divergence!(divflux, flux, space)

Compute the divergence of `flux`, storing the result in `divflux`.

Given a vector field `v`, we compute the divergence as

    [∂(Jv¹)/∂ξ¹ + ∂(Jv²)/∂ξ²] / J

where `J` is the Jacobian determinant, `vⁱ` is the `i`th contravariant component of `v`.

This is discretized at the quadrature points as

    I{[∂(I{Jv¹})/∂ξ¹ + ∂(I{Jv²})/∂ξ²] / J}

where `I{x}` is the interpolation operator applied to a field `x`.

## References
 - Taylor and Fournier (2010), equation 15
"""
function slab_divergence!(divflux, flux, space)
    # all derivatives calculated in the reference local geometry with FT precision
    FT = Spaces.undertype(space)
    D = Quadratures.differentiation_matrix(FT, space.quadrature_style)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)

    # for each element in the element stack
    Nh = length(flux)
    for h in 1:Nh
        divflux_slab = slab(divflux, h)
        flux_slab = slab(flux, h)
        local_geometry_slab = slab(space.local_geometry, h)

        ST = eltype(divflux)
        # Shared on GPU
        Jv¹ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        Jv² = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # compute flux in contravariant coordinates (v¹,v²)
            # alternatively we could do this conversion _after_ taking the derivatives
            # may have an effect on the accuracy
            # materialize if lazy
            F = flux_slab[i, j]  # materialize the flux
            Jv¹[i, j] = RecursiveApply.rmap(
                x ->
                    local_geometry.J *
                    Geometry.contravariant1(x, local_geometry),
                F,
            )
            Jv²[i, j] = RecursiveApply.rmap(
                x ->
                    local_geometry.J *
                    Geometry.contravariant2(x, local_geometry),
                F,
            )
        end
        # GPU synchronize
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # compute spectral deriv along first dimension
            ∂₁v₂ = RecursiveOperators.rmatmul1(D, Jv¹, i, j) # ∂(Jv¹)/∂ξ¹ = D[i,:]*Jv¹[:,j]
            # compute spectral deriv along second dimesion
            ∂₂v₁ = RecursiveOperators.rmatmul2(D, Jv², i, j)
            divflux_slab[i, j] =
                Contravariant3Vector(inv(local_geometry.J) ⊠ (∂₁v₂ ⊟ ∂₂v₁))
        end
    end
    return divflux
end

function slab_strong_curl!(curlvec, vec, mesh)
    # all derivatives calculated in the reference local geometry with FT precision
    FT = Meshes.undertype(mesh)
    D = Quadratures.differentiation_matrix(FT, mesh.quadrature_style)
    Nq = Quadratures.degrees_of_freedom(mesh.quadrature_style)

    # for each element in the element stack
    Nh = length(flux)
    for h in 1:Nh
        curlvec = slab(curlvec, h)
        vec_slab = slab(vec, h)
        local_geometry_slab = slab(mesh.local_geometry, h)

        ST = eltype(divflux)
        # Shared on GPU
        v₁ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        v₂ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # we convert the vec to covariant coordinates (v₁, v₂),
            # and return an orthogonal vector in contravariant coordinates v³
            v₁[i, j] = RecursiveOperators.rmap(
                x -> Geometry.covariant1(x, local_geometry),
                vec_slab,
            )
            v₂[i, j] = RecursiveOperators.rmap(
                x -> Geometry.covariant2(x, local_geometry),
                vec_slab,
            )
        end
        # GPU synchronize
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # compute spectral deriv along first dimension
            ∂₁v₂ = RecursiveOperators.rmatmul1(D, v₂, i, j)
            # compute spectral deriv along second dimension
            ∂₂v₁ = RecursiveOperators.rmatmul2(D, v₁, i, j)
            curlvec_slab[i, j] =
                Contravariant3Vector(inv(local_geometry.J) ⊠ (∂₁v₂ ⊟ ∂₂v₁))
        end
    end
    return curlvec
end

function slab_weak_curl!(curlvec, vec, mesh)
    # all derivatives calculated in the reference local geometry with FT precision
    FT = Meshes.undertype(mesh)
    D = Quadratures.differentiation_matrix(FT, mesh.quadrature_style)
    Nq = Quadratures.degrees_of_freedom(mesh.quadrature_style)

    # for each element in the element stack
    Nh = length(flux)
    for h in 1:Nh
        curlvec = slab(curlvec, h)
        vec_slab = slab(vec, h)
        local_geometry_slab = slab(mesh.local_geometry, h)

        ST = eltype(divflux)
        # Shared on GPU
        v₁ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        v₂ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # we convert the vec to covariant coordinates (v₁, v₂),
            # and return an orthogonal vector in contravariant coordinates v³
            v₁[i, j] = RecursiveOperators.rmap(
                x -> Geometry.covariant1(x, local_geometry),
                vec_slab,
            )
            v₂[i, j] = RecursiveOperators.rmap(
                x -> Geometry.covariant2(x, local_geometry),
                vec_slab,
            )
        end
        # GPU synchronize
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # compute spectral deriv along first dimension
            ∂₁v₂ = RecursiveOperators.rmatmul1(D, v₂, i, j)
            # compute spectral deriv along second dimension
            ∂₂v₁ = RecursiveOperators.rmatmul2(D, v₁, i, j)
            curlvec_slab[i, j] =
                Contravariant3Vector(inv(local_geometry.J) ⊠ (∂₁v₂ ⊟ ∂₂v₁))
        end
    end
    return curlvec
end

"""
    slab_weak_divergence!(divflux, flux, space)

Compute the right-hand side of the divergence of `flux` weakly, storing the result in `divflux`.

This computes the right-hand side of the variational problem of finding θ such that

    - ⟨ϕ, J θ⟩ = ⟨∂ϕ/∂ξ¹, J u¹⟩ + ⟨∂ϕ/∂ξ², J u²⟩

for all `ϕ` (which arises by integration by parts).

Discretely it is equivalent to

    (D₁' * W * J * u¹ + D₂' * W * J * u²)

where
 - `J` is the diagonal Jacobian matrix
 - `W` is the diagonal matrix of quadrature weights
 - `D₁` and `D₂` are the discrete derivative matrices along the first and second dimensions.
"""
function slab_weak_divergence!(divflux, flux, space)
    # all derivatives calculated in the reference local geometry with FT precision
    FT = Spaces.undertype(space)
    D = Quadratures.differentiation_matrix(FT, space.quadrature_style)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)

    # for each element in the element stack
    Nh = length(flux)
    for h in 1:Nh
        divflux_slab = slab(divflux, h)
        flux_slab = slab(flux, h)
        local_geometry_slab = slab(space.local_geometry, h)

        ST = eltype(divflux)
        # Shared on GPU
        WJv¹ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        WJv² = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # compute flux in contravariant coordinates (v¹,v²)
            # alternatively we could do this conversion _after_ taking the derivatives
            # may have an effect on the accuracy
            # materialize if lazy
            F = flux_slab[i, j]
            WJv¹[i, j] = RecursiveApply.rmap(
                x ->
                    local_geometry.WJ *
                    Geometry.contravariant1(x, local_geometry),
                F,
            )
            WJv²[i, j] = RecursiveApply.rmap(
                x ->
                    local_geometry.WJ *
                    Geometry.contravariant2(x, local_geometry),
                F,
            )
        end
        # GPU synchronize

        for i in 1:Nq, j in 1:Nq
            local_geometry = local_geometry_slab[i, j]
            # compute spectral deriv along first dimension
            Dᵀ₁WJv¹ = RecursiveApply.rmatmul1(D', WJv¹, i, j) # D'WJv¹)/∂ξ¹ = D[i,:]*Jv¹[:,j]
            # compute spectral deriv along second dimension
            Dᵀ₂WJv² = RecursiveApply.rmatmul2(D', WJv², i, j) # ∂(Jv²)/∂ξ² = D[j,:]*Jv²[i,:]
            divflux_slab[i, j] = Dᵀ₁WJv¹ ⊞ Dᵀ₂WJv²
        end
    end
    return divflux
end

function slab_weak_divergence!(
    divflux_slab::DataLayouts.DataSlab2D,
    flux_slab::DataLayouts.DataSlab2D,
    space_slab::Spaces.SpectralElementSpaceSlab,
)
    # all derivatives calculated in the reference local geometry with FT precision
    FT = Spaces.undertype(space_slab)
    D = Quadratures.differentiation_matrix(FT, space_slab.quadrature_style)
    Nq = Quadratures.degrees_of_freedom(space_slab.quadrature_style)

    # for each element in the element stack
    local_geometry_slab = space_slab.local_geometry

    ST = eltype(divflux_slab)
    # Shared on GPU
    WJv¹ = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
    WJv² = MArray{Tuple{Nq, Nq}, ST, 2, Nq * Nq}(undef)
    for i in 1:Nq, j in 1:Nq
        local_geometry = local_geometry_slab[i, j]
        # compute flux in contravariant coordinates (v¹,v²)
        # alternatively we could do this conversion _after_ taking the derivatives
        # may have an effect on the accuracy
        # materialize if lazy
        F = flux_slab[i, j]
        WJv¹[i, j] = RecursiveApply.rmap(
            x ->
                local_geometry.WJ * Geometry.contravariant1(x, local_geometry),
            F,
        )
        WJv²[i, j] = RecursiveApply.rmap(
            x ->
                local_geometry.WJ * Geometry.contravariant2(x, local_geometry),
            F,
        )
    end
    # GPU synchronize

    for i in 1:Nq, j in 1:Nq
        # compute spectral deriv along first dimension
        Dᵀ₁WJv¹ = RecursiveApply.rmatmul1(D', WJv¹, i, j) # D'WJv¹)/∂ξ¹ = D[i,:]*Jv¹[:,j]
        # compute spectral deriv along second dimension
        Dᵀ₂WJv² = RecursiveApply.rmatmul2(D', WJv², i, j) # ∂(Jv²)/∂ξ² = D[j,:]*Jv²[i,:]
        divflux_slab[i, j] = Dᵀ₁WJv¹ ⊞ Dᵀ₂WJv²
    end
    return divflux_slab
end

function slab_gradient!(∇field::Field, field::Field)
    @assert axes(∇field) === axes(field)
    Operators.slab_gradient!(
        Fields.field_values(∇field),
        Fields.field_values(field),
        axes(field),
    )
    return ∇field
end
function slab_divergence!(divflux::Field, flux::Field)
    @assert axes(divflux) === axes(flux)
    Operators.slab_divergence!(
        Fields.field_values(divflux),
        Fields.field_values(flux),
        axes(flux),
    )
    return divflux
end
function slab_weak_divergence!(divflux::Field, flux::Field)
    @assert axes(divflux) === axes(flux)
    Operators.slab_weak_divergence!(
        Fields.field_values(divflux),
        Fields.field_values(flux),
        axes(flux),
    )
    return divflux
end

function slab_gradient(field::Field)
    S = eltype(field)
    ∇S = RecursiveApply.rmaptype(T -> Cartesian12Vector{T}, S)
    Operators.slab_gradient!(similar(field, ∇S), field)
end

function slab_divergence(field::Field)
    S = eltype(field)
    divS = RecursiveApply.rmaptype(Geometry.divergence_result_type, S)
    Operators.slab_divergence!(similar(field, divS), field)
end
function slab_weak_divergence(field::Field)
    S = eltype(field)
    divS = RecursiveApply.rmaptype(Geometry.divergence_result_type, S)
    Operators.slab_weak_divergence!(similar(field, divS), field)
end

function interpolate(space_to::AbstractSpace, field_from::Field)
    field_to = similar(field_from, (space_to,), eltype(field_from))
    interpolate!(field_to, field_from)
end
function interpolate!(field_to::Field, field_from::Field)
    space_to = axes(field_to)
    space_from = axes(field_from)
    # @assert space_from.topology == space_to.topology

    M = Quadratures.interpolation_matrix(
        Float64,
        space_to.quadrature_style,
        space_from.quadrature_style,
    )
    Operators.tensor_product!(
        Fields.field_values(field_to),
        Fields.field_values(field_from),
        M,
    )
    return field_to
end

function restrict!(field_to::Field, field_from::Field)
    space_to = axes(field_to)
    space_from = axes(field_from)
    # @assert space_from.topology == space_to.topology

    M = Quadratures.interpolation_matrix(
        Float64,
        space_from.quadrature_style,
        space_to.quadrature_style,
    )
    Operators.tensor_product!(
        Fields.field_values(field_to),
        Fields.field_values(field_from),
        M',
    )
    return field_to
end

# matrix interpolate used for 2D field plots
function matrix_interpolate(
    field::Fields.SpectralElementField2D,
    Q_interp::Quadratures.Uniform{Nu},
) where {Nu}
    S = eltype(field)
    space = axes(field)
    topology = Spaces.topology(space)
    quadrature_style = Spaces.quadrature_style(space)
    mesh = topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2
    interp_data =
        DataLayouts.IH1JH2{S, Nu}(Matrix{S}(undef, (Nu * n1, Nu * n2)))
    M = Quadratures.interpolation_matrix(Float64, Q_interp, quadrature_style)
    Operators.tensor_product!(interp_data, Fields.field_values(field), M)
    return parent(interp_data)
end

function matrix_interpolate(
    field::Fields.ExtrudedFiniteDifferenceField,
    Q_interp::Quadratures.Uniform{Nu},
) where {Nu}
    S = eltype(field)
    space = axes(field)
    quadrature_style = Spaces.quadrature_style(space)
    nl = Spaces.nlevels(space)
    mesh = Spaces.topology(space).mesh
    n1 = mesh.n1
    interp_data = DataLayouts.IV1JH2{S, Nu}(Matrix{S}(undef, (nl, Nu * n1)))
    M = Quadratures.interpolation_matrix(Float64, Q_interp, quadrature_style)
    Operators.tensor_product!(interp_data, Fields.field_values(field), M)
    return parent(interp_data)
end

matrix_interpolate(field::Field, Nu::Integer) =
    matrix_interpolate(field, Quadratures.Uniform{Nu}())


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 774bb4f2-963e-11ec-30ff-e3394f52858b
using HypertextLiteral, Parameters, PlutoDevMacros, Deno_jll

# ╔═╡ d12da0cd-80fe-4c09-907b-d20a81c07c1c
@only_in_nb using PlutoUI

# ╔═╡ 320477a5-77ad-40b3-bccd-37a06c06c22e
@only_in_nb TableOfContents()

# ╔═╡ b1d3d2ca-aa67-47de-a2a8-6119f3ca4e46
md"# Ingredients"

# ╔═╡ bad43269-199b-4176-8548-8a963268cb46
@plutoinclude "./js_base.jl" "all"

# ╔═╡ 09e4ee53-b8b7-4de0-9439-d1bc7bc22578
md"# Javascript Snippet macro"

# ╔═╡ 4d193bff-9ed2-42ba-9cef-4d8c9abc8332
esc_if_needed(x::String) = x

# ╔═╡ c9458743-0471-430e-b7a5-e26a2ccb5811
esc_if_needed(x) = esc(x)

# ╔═╡ 279012f1-b67a-4fac-a169-34dbe47e4638
begin
	struct RenderWithoutScriptTags
		x
	end

	function Base.show(io::IO, m::MIME"text/javascript", r::RenderWithoutScriptTags)
		full_result = repr(MIME"text/html"(), r.x; context=io)
		write(io,
			"((ctx, x_scale, y_scale) => {" * full_result[1+length("<script>"):end-length("</script>")] * "})"
		)
	end
end

# ╔═╡ 227f7480-a39d-4692-9874-ff9328e4080a
begin
macro js(str_expr::Expr)
	QuoteNode(str_expr)
	quote
		result = @htl( $(Expr(:string, "<script>", esc_if_needed.(str_expr.args)..., "</script>")) )
		RenderWithoutScriptTags(result)
	end
end
macro js(str_expr::String)
	QuoteNode(str_expr)

	quote
		result = @htl( $(Expr(:string, "<script>", esc_if_needed(str_expr), "</script>")) )
		RenderWithoutScriptTags(result)
	end
end
end

# ╔═╡ e3156da9-8603-41f9-ad41-03eaf79c4540
md"# D3"

# ╔═╡ b5110c15-0b46-4450-b8d2-9a83f7aeef54
md"## D3 Component Type"

# ╔═╡ c06a51f5-224a-4ffa-87a7-d5f39a1a4bb5
abstract type D3Component end

# ╔═╡ 37f3557c-0d82-464d-9f36-847eceebaae6
Base.extrema(component::D3Component) =
	error("Extrema not implemented for type $(typeof(component))")

# ╔═╡ d8c59da2-5a98-4321-89f2-4766d9ad14c6
md"## D3 Attributes"

# ╔═╡ c80d5063-aac3-4ee7-994c-eb6394b225eb
@with_kw struct D3Attr
  attr     :: NamedTuple = (;)
  style    :: NamedTuple = (;)
  duration :: Int = 200
  events   :: Vector{<:AbstractString} = String[]
end

# ╔═╡ 813d7b17-a806-4d9e-95ed-e095a43a8291
to_named_tuple(attr::D3Attr) = (attr=attr.attr, style=attr.style, duration=attr.duration, events=attr.events)

# ╔═╡ 1eaa1bb9-7893-4181-973c-feda25677ef3
md"## Curve Types"

# ╔═╡ 25710a40-fe6a-4f1d-8851-00ddc810d8f4
@enum Curve Cardinal Natural CatmullRom MonotoneX MonotoneY Basis BasisClosed Linear LinearClosed

# ╔═╡ 5e2b8127-3385-406a-8de1-4acb74ff3126
md"## Shape Types"

# ╔═╡ e19d9ac1-b314-4901-a448-f58230011b05
@enum ShapeType Square Circle Triangle

# ╔═╡ 9763e4b9-515d-4576-9f76-1b337365d6a6
Base.show(io::IO,  m::MIME"text/javascript", shapeType::ShapeType) = Base.show(io, m, HypertextLiteral.JavaScript("d3.symbol" * (string(shapeType))))

# ╔═╡ 743f182a-46da-4732-b1e4-5a237912848f
Base.show(io::IO,  m::MIME"text/javascript", curve::Curve) = Base.show(io, m, HypertextLiteral.JavaScript("curve" * string(curve)))

# ╔═╡ 8803b1d9-a29c-4b4c-b5bc-93a78561a70f
begin
	struct PublishToJS
	x
	end
	function Base.show(io::IO, m::MIME"text/javascript", ptj::PublishToJS)
		Base.show(io, m, Main.PlutoRunner.publish_to_js(ptj.x))
	end
	better_publish_to_js(x) = PublishToJS(x)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Deno_jll = "04572ae6-984a-583e-9378-9577a1c2574d"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
PlutoDevMacros = "a0499f29-c39b-4c5c-807c-88074221b949"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Deno_jll = "~1.20.4"
HypertextLiteral = "~0.9.3"
Parameters = "~0.12.3"
PlutoDevMacros = "~0.4.5"
PlutoUI = "~0.7.38"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Deno_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "970da1e64a94f13b51c81691c376a1d5a83a0b3c"
uuid = "04572ae6-984a-583e-9378-9577a1c2574d"
version = "1.20.4+0"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoDevMacros]]
deps = ["MacroTools", "Requires"]
git-tree-sha1 = "994167def8f46d3be21783a76705228430e29632"
uuid = "a0499f29-c39b-4c5c-807c-88074221b949"
version = "0.4.5"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═774bb4f2-963e-11ec-30ff-e3394f52858b
# ╟─d12da0cd-80fe-4c09-907b-d20a81c07c1c
# ╟─320477a5-77ad-40b3-bccd-37a06c06c22e
# ╟─b1d3d2ca-aa67-47de-a2a8-6119f3ca4e46
# ╠═bad43269-199b-4176-8548-8a963268cb46
# ╟─09e4ee53-b8b7-4de0-9439-d1bc7bc22578
# ╠═4d193bff-9ed2-42ba-9cef-4d8c9abc8332
# ╠═c9458743-0471-430e-b7a5-e26a2ccb5811
# ╠═279012f1-b67a-4fac-a169-34dbe47e4638
# ╠═227f7480-a39d-4692-9874-ff9328e4080a
# ╠═8803b1d9-a29c-4b4c-b5bc-93a78561a70f
# ╟─e3156da9-8603-41f9-ad41-03eaf79c4540
# ╟─b5110c15-0b46-4450-b8d2-9a83f7aeef54
# ╠═c06a51f5-224a-4ffa-87a7-d5f39a1a4bb5
# ╠═37f3557c-0d82-464d-9f36-847eceebaae6
# ╟─d8c59da2-5a98-4321-89f2-4766d9ad14c6
# ╠═c80d5063-aac3-4ee7-994c-eb6394b225eb
# ╠═813d7b17-a806-4d9e-95ed-e095a43a8291
# ╟─1eaa1bb9-7893-4181-973c-feda25677ef3
# ╠═25710a40-fe6a-4f1d-8851-00ddc810d8f4
# ╠═743f182a-46da-4732-b1e4-5a237912848f
# ╟─5e2b8127-3385-406a-8de1-4acb74ff3126
# ╠═e19d9ac1-b314-4901-a448-f58230011b05
# ╠═9763e4b9-515d-4576-9f76-1b337365d6a6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
using eFEMpart
include("Stokes_example_param.jl")

function stokes_example(p::Stokes_Para)
  mesh = squareMeshFluid([p.xmin,p.xmax,p.ymin,p.ymax],p.Nintervals)

  param = FluidParam(p.μ)
  OperatorType = :Stokes2D

  dNodes = Dirichlet(:left,:bottom,:top)
  nNodes = Neumann(:right)
  Nodes = [dNodes,nNodes]

  dUBCf = Dirichlet((x,y) -> (x==p.xmin ? -p.lBCval*(p.ymax-y)*(p.ymin-y) : 0.0))
  dVBCf = Dirichlet((x,y) -> 0.0)
  bcfun = [dUBCf,dVBCf]

  prob = Problem(mesh,Nodes,bcfun,OperatorType)
  sol = solve(prob,mesh,param)

  # save solution
  vtkname = Path(p.solFolder,p.solFile)
  sd = ScalarData(sol.p)
  sn = ScalarNames("pressure")
  vd = VectorData([sol.u,sol.v])
  vn = VectorNames("velocity")
  vtksave(mesh,sd,sn,vd,vn,vtkname)

  nothing
end

function main()
  p = Stokes_Para()
  stokes_example(p)
end

main()

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
@testset "Limited" begin
    @testset "Limited Array" begin
        L1 = Limited([1, 2, 3]; limit = 5)

        @test n_dropped(L1) == 0

        for i in 4:8
            push!(L1, i)
        end

        @test n_dropped(L1) == 3

        L1[1] = 0

        @test L1[1] == 0
    end

    @testset "Limited Dict" begin
        L1 = Limited(Dict("a" => 1, "b" => 2); limit = 1)

        @test n_dropped(L1) == 1
        @test length(pairs(L1)) == 1
        @test length(collect(L1)) == 1

        L1["c"] = 3

        @test haskey(L1, "c") == false
        @test_throws KeyError L1["c"]
        @test n_dropped(L1) == 2
    end
end

@testset "StaticAttrs" begin
    A = StaticAttrs()
    @test length(A) == 0

    A1 = StaticAttrs((x = 1, y = -2.0, z = "abcdefg"); value_length_limit = 3)
    @info "Show StaticAttrs" A1
    @test A1["x"] == 1
    @test A1["z"] == "abc"

    A2 = StaticAttrs(x = 1, y = -2.0, z = "0")
    @test keys(A2) == (:x, :y, :z)
    @test haskey(A2, :z)
    @test A2[(:x, :y)] == StaticAttrs(x = 1, y = -2.0)
    @test A2[:z] == "0"
    @test length(collect(A2)) == 3
    @test collect(pairs(A2)) == [:x => 1, :y => -2.0, :z => "0"]

    A3 = StaticAttrs("test.x" => 1, "test.y" => -2.0, "test.z" => "0")
    @test A3["test.z"] == "0"

    @test_throws ArgumentError StaticAttrs("abc" => nothing)  # invalid value type

    # always no elements dropped
    @test n_dropped(A) == 0
    @test n_dropped(A1) == 0
    @test n_dropped(A2) == 0

    S = StaticAttrs(ab = 102, bc = 101, ac = 101)
    @test sort(S) == StaticAttrs(ab = 102, ac = 101, bc = 101)
end

@testset "DynamicAttrs" begin
    d = DynamicAttrs("a" => 1; count_limit = 4, value_length_limit = 5)
    @test haskey(d, "a")
    d["b"] = "bbbbb"
    d["c"] = "ccccc" * "extra"
    d["d"] = ["ddddd" * "extra", "dd"]
    d["e"] = "???"

    @test length(d) == 4
    @test length(collect(d)) == 4
    @test length(pairs(d)) == 4

    @test n_dropped(d) == 1
    @test d["b"] == "bbbbb"
    @test d["c"] == "ccccc"
    @test d["d"] == ["ddddd", "dd"]

    d["b"] = "bbbbbbbbbb"
    @test d["b"] == "bbbbb"

    # count_limit exceeded
    d["x"] = "abcdefg"
    @test_throws KeyError d["x"]
end

@testset "context" begin
    @test isnothing(get(current_context(), :say, nothing))

    with_context(; say = "foo") do
        @test current_context()[:say] == "foo"
        with_context(; say = "bar") do
            @test current_context()[:say] == "bar"
        end
        @test current_context()[:say] == "foo"
    end

    @test isnothing(get(current_context(), :say, nothing))

    with_context(; a = "foo") do
        with_context(; b = "bar") do
            @test current_context()[:a] == "foo"
            @test current_context()[:b] == "bar"
        end
    end
end


text:
using ReachabilityAnalysis
import DataStructures
using DataStructures: OrderedDict
using JuMP, Ipopt

include("models.jl")

# print the step size used in each iteration
VERBOSE = false

# scaling factor for time unit (base is in seconds)
unit_scale = 1e3  # milliseconds

# methods
algs_homog = [
        (FirstOrderddt(), "d-dt"),
        (FirstOrderZonotope(), "Zonotope"),
        (CorrectionHull(order=4), "Correction-hull"),
        (FirstOrder(), "First-order"),
        (ForwardBackward(solver=Ipopt.Optimizer()), "Forward-backward"),
        (Forward(), "Forward")
       ]
# excludes the "d/dt" method
algs_heterog = [
        (FirstOrderZonotope(), "Zonotope"),
        (CorrectionHull(order=4), "Correction-hull"),
        (FirstOrder(), "First-order"),
        (ForwardBackward(solver=Ipopt.Optimizer()), "Forward-backward"),
        (Forward(), "Forward")
       ]

function get_time(collection, m)
    if m != 1
        # m is the number of runs per experiment
        sort!(collection)
        part = ceil(Int, m / 20)
        collection = collection[part:m-part]
    end
    time = sum(collection) / length(collection) * unit_scale
end

# run experiment
function experiment(P, δs, model_name, algs; m::Int=1)
    d = ones(dim(P.x0))  # direction for comparison
    δ2name2results = OrderedDict{Float64, OrderedDict}()
    for δ in δs
        local name2results = OrderedDict{String, Tuple}()
        for (alg, name) in algs
            collection = Vector{Float64}(undef, m)
            local Ω0
            @inbounds for i in 1:m
                res = @timed discretize(P, δ, alg)
                Ω0 = res.value.x0
                collection[i] = res.time
            end
            local time1 = get_time(collection, m)

            collection = Vector{Float64}(undef, m)
            local support_value
            @inbounds for i in 1:m
                res = @timed ρ(d, Ω0)
                support_value = res.value
                collection[i] = res.time
            end
            local time2 = get_time(collection, m)
            local tuple = (Ω0, support_value, time1 + time2)
            name2results[name] = tuple
        end
        δ2name2results[δ] = name2results
        VERBOSE && println("δ = $δ")
    end

    # write results to table
    for (idx, filename) in [(2, "delta_support_$(model_name).dat"),
                            (3, "delta_time_$(model_name).dat")]
        open(filename, "w") do f
            # header
            write(f, "delta")
            for (_, name) in algs
                write(f, "\t$name")
            end
            write(f, "\n")

            # data
            for (δ, name2results) in δ2name2results
                write(f, "$δ")
                for (i, (name, tuple)) in enumerate(name2results)
                    @assert name == algs[i][2] "$name != $(algs[i][2])"
                    data = tuple[idx]
                    write(f, "\t$data")
                end
                write(f, "\n")
            end
        end
    end

    return δ2name2results
end

# delta steps
δs_oscillator = vcat(
          range(0.01, stop=0.06, step=0.01)
         )
δs_freedom = vcat(
          range(0.000005, stop=0.0001, step=0.000005),
          range(0.0002, stop=0.001, step=0.0001),
          range(0.002, stop=0.01, step=0.001)
         )
δs_iss = vcat(
          [0.00002],  # precompilation
          range(0.00002, stop=0.0001, step=0.00001),
          range(0.0002, stop=0.001, step=0.0001),
          range(0.002, stop=0.02, step=0.001)
         )

# run experiments
println("Running oscillator model ...")
experiment(oscillator(), δs_oscillator, "oscillator", algs_homog; m=40)
println("Running sdof model ...")
experiment(freedom(), δs_freedom, "freedom", algs_homog; m=40)
println("Running ISS model ...")
experiment(iss(), δs_iss, "iss", algs_heterog)

nothing


text:
using AbstractGPs
using GPARs
using Random
using Test

using GPARs: posterior

@testset "GPARs.jl" begin
    include("gpar.jl")
end


text:
using NearestNeighbors: KDTree
@testset "in 2d" begin
    h = 0.3
    Δx = 0.03
    Δy = 0.05
    x = mapreduce(identity, hcat, [[x; y] for x in 0 : Δx : 1 for y in 0 : Δy : 2])
    diams = 2 * norm([Δx, Δy]) * ones(size(x, 2))
    domain_vector = CompressingSolvers.array2domains(x)
    @test all(CompressingSolvers.iselementary.(domain_vector))
    # Construct a new "parent domain" that contains all children domains
    parent_domain = CompressingSolvers.Domain(domain_vector, length(domain_vector) + 1)
    # creating the hierarchy from the parent domain
    hierarchy = CompressingSolvers.create_hierarchy((CompressingSolvers.children(parent_domain)), h, diams, KDTree)
    # returns an array of arrays of domains
    all_domains = CompressingSolvers.gather_hierarchy(hierarchy)
    for k = 1 : length(all_domains)
        sorted_truth = sort(domain_vector, by=CompressingSolvers.id)
        sorted_recovered = sort(CompressingSolvers.gather_descendants(all_domains[k]), by=CompressingSolvers.id)

        # make sure the weights stay the same
        @test sum(CompressingSolvers.weight.(all_domains[k])) == sum(CompressingSolvers.weight.(domain_vector))

        # Test if the partition, on each level, produces the same set of elements
        @test sort(domain_vector, by=CompressingSolvers.id) == sort(CompressingSolvers.gather_descendants(all_domains[k]), by=CompressingSolvers.id)
    end
    CompressingSolvers.plot_domains(all_domains[1]; xlims=(0.0, 2.0), ylims=(0.0,2.0))
end

@testset "in 2d, inhomogeneous" begin
    h = 0.3
    Δratio = 3
    Δx = 0.03
    Δy = 0.05
    x1 = mapreduce(identity, hcat, [[x; y] for x in 0 : Δx : 1 for y in 0 : Δy : 2])
    x2 = mapreduce(identity, hcat, [[x; y] for x in (1 + Δx * Δratio)  : (Δx * Δratio) : 2 for y in (Δy * Δratio): (Δy * Δratio) : 1])
    diams = vcat(2 * norm([Δx, Δy]) * ones(size(x2, 2)) * Δratio, 2 * norm([Δx, Δy]) * ones(size(x1, 2)))
    x = hcat(x2, x1)
    domain_vector = CompressingSolvers.array2domains(x, vcat(Δratio^2 * ones(size(x2, 2)), ones(size(x1, 2))))
    @test all(CompressingSolvers.iselementary.(domain_vector))
    # Construct a new "parent domain" that contains all children domains
    parent_domain = CompressingSolvers.Domain(domain_vector, length(domain_vector) + 1)
    # creating the hierarchy from the parent domain
    hierarchy = CompressingSolvers.create_hierarchy((CompressingSolvers.children(parent_domain)), h, diams, KDTree)
    # returns an array of arrays of domains
    all_domains = CompressingSolvers.gather_hierarchy(hierarchy)

    # Test whether there are no "spurious scales"
    @test length(unique(length.(all_domains))) == length(all_domains)
    for k = 1 : length(all_domains)
        sorted_truth = sort(domain_vector, by=CompressingSolvers.id)
        sorted_recovered = sort(CompressingSolvers.gather_descendants(all_domains[k]), by=CompressingSolvers.id)

        if k < 3
            # make sure the weights stay the same
            @test sum(CompressingSolvers.weight.(all_domains[k])) == sum(CompressingSolvers.weight.(domain_vector))

            # Test if the partition, on each level, produces the same set of elements
            @test sort(domain_vector, by=CompressingSolvers.id) == sort(CompressingSolvers.gather_descendants(all_domains[k]), by=CompressingSolvers.id)
        end
    end
    CompressingSolvers.plot_domains(all_domains[1]; xlims=(0.0, 2.0), ylims=(0.0,2.0))
end

@testset "in 3d" begin
    Δx = 0.06
    Δy = 0.07
    Δz = 0.08
    x = mapreduce(identity, hcat, [[x; y; z] for x in 0 : Δx : 1 for y in 0 : Δy : 2 for z in -0.3 : Δz : 0.2])
    diams = 2 * norm([Δx, Δy, Δz]) * ones(size(x, 2))
    h = 0.5
    domain_vector = CompressingSolvers.array2domains(x)
    @test all(CompressingSolvers.iselementary.(domain_vector))
    # Construct a new "parent domain" that contains all children domains
    parent_domain = CompressingSolvers.Domain(domain_vector, length(domain_vector) + 1)
    # creating the hierarchy from the parent domain
    hierarchy = CompressingSolvers.create_hierarchy((CompressingSolvers.children(parent_domain)), h, diams, KDTree)
    # returns an array of arrays of domains
    all_domains = CompressingSolvers.gather_hierarchy(hierarchy)
    # Test whether there are no "spurious scales"
    @test length(unique(length.(all_domains))) == length(all_domains)
    for k = 1 : length(all_domains)
        sorted_truth = sort(domain_vector, by=CompressingSolvers.id)
        sorted_recovered = sort(CompressingSolvers.gather_descendants(all_domains[k]), by=CompressingSolvers.id)

        # Test if the partition, on each level, produces the same set of elements
        @test sort(domain_vector, by=CompressingSolvers.id) == sort(CompressingSolvers.gather_descendants(all_domains[k]), by=CompressingSolvers.id)
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    kekulize,
    trivialhydrogens, allhydrogens,
    removehydrogens, addhydrogens,
    largestcomponentnodes, extractlargestcomponent,
    protonateacids, deprotonateoniums,
    depolarize, polarize, totriplebond, toallenelike

# TODO: large conjugated system
# TODO: salts and waters should detected by functional group analysis
# TODO: Phosphate, diphosphate, sulfate, nitrate, acetate,
# maleate, fumarate, succinate, citrate, tartrate, oxalate,
# mesylate, tosylate, besylate,
# benzoate, gluconate

"""
    kekulize(mol::SMILES) -> SMILES

Return the molecule with SMILES aromatic bonds kekulized.

SMILES allows aromatic atoms in small letters - b, c, n, o, p, s, [as] and [se]. Once these are stored in `SmilesAtom.isaromatic` field, then `kekulize` will place double bonds to satisfy valences.

Kekulization should be applied for molecules parsed from SMILES. Otherwise, some bond valence and implicit hydrogen properties would be wrong.
"""
function kekulize(mol::SMILES)
    nodes = Set{Int}()
    for i in 1:nodecount(mol)
        nodeattr(mol, i).isaromatic === nothing && continue
        nodeattr(mol, i).isaromatic || continue
        if atomsymbol(mol)[i] === :C
            push!(nodes, i)
        elseif atomsymbol(mol)[i] in (:N, :P, :As)
            if (degree(mol, i) == 2 ||
                (degree(mol, i) == 3 && nodeattr(mol, i).charge == 1))
                push!(nodes, i)
            end
        end
    end
    subg = nodesubgraph(mol, nodes)
    a, b = twocoloring(subg)
    adjmap = Dict{Int,Set{Int}}(
        i => adjacencies(subg, i) for i in nodeset(subg))
    mapping = maxcardmap(a, b, adjmap)
    newmol = graphmol(mol)
    for (u, v) in mapping
        e = findedgekey(mol, u, v)
        setedgeattr!(newmol, e, setorder(edgeattr(newmol, e), 2))
    end
    return newmol
end

"""
    trivialhydrogens(mol::GraphMol) -> Set{Int}

Return a set of trivial hydrogen nodes.

"Trivial" means not charged, no unpaired electron, no specific mass, non-stereospecific and attached to organic heavy atoms. Note that stereochemstry of SDFile molecules will not be considered until `setstereocenter!` is applied.
"""
function trivialhydrogens(mol::GraphMol)
    hs = Set{Int}()
    organic_heavy = (
        :B, :C, :N, :O, :F, :Si, :P, :S, :Cl, :As, :Se, :Br, :I)
    for (i, atom) in enumerate(nodeattrs(mol))
        atom.symbol == :H || continue
        atom.charge == 0 || continue
        atom.multiplicity == 1 || continue
        atom.mass === nothing || continue
        @assert atom.stereo === :unspecified
        degree(mol, i) == 1 || continue
        (inc, adj) = iterate(neighbors(mol, i))[1]
        nodeattr(mol, adj).symbol in organic_heavy || continue
        edgeattr(mol, inc).order == 1 || continue
        nodeattr(mol, adj).stereo === :unspecified || continue
        if nodeattr(mol, adj) isa SmilesAtom
            @assert !nodeattr(mol, adj).isaromatic
        end
        push!(hs, i)
    end
    return hs
end

"""
    allhydrogens(mol::GraphMol) -> Set{Int}

Return a set of all hydrogen nodes.
"""
function allhydrogens(mol::GraphMol)
    hs = Set{Int}()
    for (i, atom) in enumerate(nodeattrs(mol))
        atom.symbol == :H || continue
        push!(hs, i)
    end
    return hs
end

"""
    removehydrogens(mol::GraphMol) -> GraphMol

Return the molecule with hydrogen nodes removed.

If option `all` is set to true (default), all hydrogens will be removed, otherwise only trivial hydrogens will be removed (see [`trivialhydrogens`](@ref)).
"""
function removehydrogens(mol::GraphMol; all=true)
    hydrogens = all ? allhydrogens : trivialhydrogens
    ns = setdiff(nodeset(mol), hydrogens(mol))
    return graphmol(nodesubgraph(mol, ns))
end

"""
    addhydrogens(mol::GraphMol) -> GraphMol

Return the molecule with all hydrogen nodes explicitly attached.
"""
function addhydrogens(mol::GraphMol{A,B}) where {A<:Atom,B<:Bond}
    implicithcount_ = implicithcount(mol)
    newmol = graphmol(mol)
    for i in 1:nodecount(mol)
        for j in 1:implicithcount_[i]
            addedge!(newmol, i, addnode!(newmol, A(:H)), B())
        end
    end
    return newmol
end

"""
    largestcomponentnodes(mol::GraphMol) -> Set{Int}

Return a set of nodes in the largest connected component.
"""
largestcomponentnodes(mol::GraphMol
    ) = sortstablemax(connected_components(mol), by=length, init=Set{Int}())

"""
    extractlargestcomponent(mol::GraphMol) -> GraphMol

Return the largest connected component of the molecular graph.

This should be useful when you want to remove salt and water molecules from the molecular graph simply. On the other hand, this can remove important components from the mixture so carefully apply this preprocess method.
"""
extractlargestcomponent(mol::GraphMol
    ) = graphmol(nodesubgraph(mol, largestcomponentnodes(mol)))

"""
    protonateacids(mol::GraphMol) -> GraphMol

Return the molecule with its oxo(thio) anion groups protonated.
"""
function protonateacids(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    connectivity_ = connectivity(mol)
    newmol = graphmol(mol)
    for o in findall(charge_ .== -1)
        atomsymbol_[o] in (:O, :S) || continue
        @assert connectivity_[o] == 1
        nbr = pop!(adjacencies(mol, o))
        charge_[nbr] == 1 && continue  # polarized double bond
        setnodeattr!(newmol, o, setcharge(nodeattr(newmol, o), 0))
    end
    return newmol
end

"""
    deprotonateoniums(mol::GraphMol) -> GraphMol

Return the molecule with its onium groups are deprotonated.
"""
function deprotonateoniums(mol::GraphMol)
    hcount_ = hcount(mol)
    charge_ = charge(mol)
    newmol = graphmol(mol)
    for o in findall(charge_ .== 1)
        hcount_[o] > 0 || continue
        setnodeattr!(newmol, o, setcharge(nodeattr(newmol, o), 0))
    end
    return newmol
end

"""
    depolarize(mol::GraphMol; negative=:O, positive=[:C, :P]) -> GraphMol

Return the molecule with its dipole double bonds are depolarized.
"""
function depolarize(mol::GraphMol; negative=:O, positive=[:C, :P])
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    connectivity_ = connectivity(mol)
    isaromatic_ = isaromatic(mol)
    newmol = graphmol(mol)
    for o in findall(charge_ .== -1)
        atomsymbol_[o] === negative || continue
        @assert connectivity_[o] == 1
        (inc, adj) = iterate(neighbors(mol, o))[1]
        atomsymbol_[adj] in positive || continue
        charge_[adj] == 1 || continue
        isaromatic_[adj] && continue
        setnodeattr!(newmol, o, setcharge(nodeattr(newmol, o), 0))
        setnodeattr!(newmol, adj, setcharge(nodeattr(newmol, adj), 0))
        setedgeattr!(newmol, inc, setorder(edgeattr(newmol, inc), 2))
    end
    return newmol
end

"""
    polarize(mol::GraphMol; negative=:O, positive=[:N, :S]) -> GraphMol

Return the molecule with its dipole double bonds are polarized.
"""
function polarize(mol::GraphMol, negative=:O, positive=[:N, :S])
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    bondorder_ = bondorder(mol)
    connectivity_ = connectivity(mol)
    newmol = graphmol(mol)
    for o in findall(connectivity_ .== 1)
        atomsymbol_[o] === negative || continue
        charge_[o] == 0 || continue
        (inc, adj) = iterate(neighbors(mol, o))[1]
        atomsymbol_[adj] in positive || continue
        charge_[adj] == 0 || continue
        bondorder_[inc] == 2 || continue
        setnodeattr!(newmol, o, setcharge(nodeattr(newmol, o), -1))
        setnodeattr!(newmol, adj, setcharge(nodeattr(newmol, adj), 1))
        setedgeattr!(newmol, inc, setorder(edgeattr(newmol, inc), 1))
    end
    return newmol
end

function find13dipoles(mol::GraphMol)
    charge_ = charge(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    triads = Tuple{Int,Int,Int}[]
    for c in findall((pie_ .== 2) .* (charge_ .== 1))
        negs = Int[]
        notnegs = Int[]
        for (inc, adj) in neighbors(mol, c)
            push!(charge_[adj] == -1 ? negs : notnegs, adj)
        end
        (length(negs) == 1 && length(notnegs) == 1) || continue
        push!(triads, (negs[1], c, notnegs[1]))
    end
    return triads
end

"""
    totriplebond(mol::GraphMol) -> GraphMol

Return the molecule with its 1,3-dipole groups described in triple bond resonance structure (ex. Diazo group C=[N+]=[N-] -> [C-][N+]#N).
"""
function totriplebond(mol::GraphMol)
    newmol = graphmol(mol)
    for (first, center, third) in find13dipoles(mol)
        fe = findedgekey(mol, first, center)
        te = findedgekey(mol, center, third)
        edgeattr(mol, fe).order == 2 || continue
        setnodeattr!(newmol, first, setcharge(nodeattr(mol, first), 0))
        setnodeattr!(newmol, third, setcharge(nodeattr(mol, third), -1))
        setedgeattr!(newmol, fe, setorder(edgeattr(mol, fe), 3))
        setedgeattr!(newmol, te, setorder(edgeattr(mol, te), 1))
    end
    return newmol
end

"""
    toallenelike(mol::GraphMol) -> GraphMol

Return the molecule with its 1,3-dipole groups described in allene-like resonance structure (ex. Diazo group [C-][N+]#N -> C=[N+]=[N-]).
"""
function toallenelike(mol::GraphMol)
    newmol = graphmol(mol)
    for (first, center, third) in find13dipoles(mol)
        fe = findedgekey(mol, first, center)
        te = findedgekey(mol, center, third)
        edgeattr(mol, fe).order == 1 || continue
        setnodeattr!(newmol, first, setcharge(nodeattr(mol, first), 0))
        setnodeattr!(newmol, third, setcharge(nodeattr(mol, third), -1))
        setedgeattr!(newmol, fe, setorder(edgeattr(mol, fe), 2))
        setedgeattr!(newmol, te, setorder(edgeattr(mol, te), 2))
    end
    return newmol
end


text:
function validate(mesh::PolytopeVertexMesh)
    # non-negative measure,
    # no overlap of edges
    # convexity of linear meshes
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using IcanProj
using JFSM
using PyPlot

function run_new(filename, year, month, day, hour, dt)

    rainf, snowf, tair, iswr, ilwr, pres, rhum, wind = read_mtclim_fsm(filename)

    metdata = [year month day hour iswr ilwr snowf/dt rainf/dt tair+273.15 rhum wind pres*1000]

    am = 0
    cm = 0
    dm = 0
    em = 0
    hm = 0

    md = FsmType(am, cm, dm, em, hm, dt)

    hs = run_fsm(md, metdata)

    return hs

end

filename = "../data/metdata_71.16148_25.54974"

dt = 3*3600

time = [DateTime(2002,1,1,6)+i*Dates.Hour(3) for i in 0:(length(rainf)-1)]

year = Dates.year.(time)

month = Dates.month.(time)

day = Dates.day.(time)

hour = Dates.hour.(time)

hs = run_new(filename, year, month, day, hour, dt)


text:
"""
```
distance(x,y)
```

Functions to compute Euclidean distances between 2 3-dimensional vectors.

"""
@inline distance(x, y) = sqrt((y[1] - x[1])^2 + (y[2] - x[2])^2 + (y[3] - x[3])^2)

"""

```
distance(x,y,sides)
```

Functions to compute Euclidean distances between two 3-dimensional vectors,
with periodic boundary conditions given by `sides`.

"""
@inline function distance(x, y, sides)
    d = 0.0
    for i = 1:3
        dx = (x[i] - y[i]) % sides[i]
        if dx > sides[i] / 2
            dx = dx - sides[i]
        elseif dx < -sides[i] / 2
            dx = dx + sides[i]
        end
        d += dx^2
    end
    sqrt(d)
end


text:
ti = time()

# Continuous and discrete
include("lyapunov_exponents.jl")
include("gali_tests.jl")
include("partially_predictable_tests.jl")

# Continuous
include("poincare_tests.jl")

# Discrete
include("orbitdiagram_tests.jl")
include("periodicity_tests.jl")

# Numeric
include("entropy_dimension.jl")
include("nlts_tests.jl")
include("period_tests.jl")

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits = 33), " seconds or ",
round(ti/60, digits = 3), " minutes")


text:
using CompScienceMeshes
using BEAST

using Test

o, x, y, z = euclidianbasis(3)
tet = simplex(x,y,z,o)

nbd1 = neighborhood(tet, [0,1,1]/3)
nbd2 = neighborhood(tet, [1,0,1]/3)
nbd3 = neighborhood(tet, [1,1,0]/3)
nbd4 = neighborhood(tet, [1,1,1]/3)

rs = BEAST.NDLCDRefSpace{Float64}()
fcs = BEAST.faces(tet)
@test dot(rs(nbd1)[1].value, normal(fcs[1])) > 0
@test dot(rs(nbd2)[2].value, normal(fcs[2])) > 0
@test dot(rs(nbd3)[3].value, normal(fcs[3])) > 0
@test dot(rs(nbd4)[4].value, normal(fcs[4])) > 0

ctr = cartesian(center(tet))
@test dot(cartesian(nbd1)-ctr, normal(fcs[1])) > 0
@test dot(cartesian(nbd2)-ctr, normal(fcs[2])) > 0
@test dot(cartesian(nbd3)-ctr, normal(fcs[3])) > 0
@test dot(cartesian(nbd4)-ctr, normal(fcs[4])) > 0


text:
module QIntegrations

using Reexport
@reexport using QuadGK
#@reexport using HCubature
#@reexport using Cuba

export integrate_1d

#include("scipy_integration.jl")
#include("cuba_integration.jl")
#include("hcubature_integration.jl")
include("quadgk_integration.jl")
include("cpvagk.jl")

"""
    integrate_1d(f, a, b; args=(), rtol=1e-6, atol=1e-8)

Calculate numerical integration of function `f` from `a` to `b`.
"""
function integrate_1d(f, a, b; args=(), rtol=1e-6, atol=1e-8, kwargs...)
    integrand, lower_limit, upper_limit = integral_limit_transfromation(f, a, b, args)
    quadgk(integrand, lower_limit, upper_limit, rtol=rtol, atol=atol; kwargs...)
end

end # module


text:
module TimeoutRequest

using ..ConnectionPool
import ..@debug, ..DEBUG_LEVEL
import ..Streams: Stream

struct ReadTimeoutError <:Exception
    readtimeout::Int
end

function Base.showerror(io::IO, e::ReadTimeoutError)
    print(io, "ReadTimeoutError: Connection closed after $(e.readtimeout) seconds")
end

export timeoutlayer

"""
    timeoutlayer(stream) -> HTTP.Response

Close `IO` if no data has been received for `timeout` seconds.
"""
function timeoutlayer(handler)
    return function(stream::Stream; readtimeout::Int=0, kw...)
        if readtimeout <= 0
            # skip
            return handler(stream; kw...)
        end
        io = stream.stream
        wait_for_timeout = Ref{Bool}(true)
        timedout = Ref{Bool}(false)

        @async while wait_for_timeout[]
            if isreadable(io) && inactiveseconds(io) > readtimeout
                timedout[] = true
                close(io)
                @debug 1 "💥  Read inactive > $(readtimeout)s: $io"
                break
            end
            sleep(readtimeout / 10)
        end

        try
            return handler(stream; kw...)
        catch e
            if timedout[]
            throw(ReadTimeoutError(readtimeout))
            end
            rethrow(e)
        finally
            wait_for_timeout[] = false
        end
    end
end

end # module TimeoutRequest


text:
"""
    mutable struct Structure_Ti  <: Str

This contains information for model constructing of Ti and its compounds.

# Supported list:
- structure_tio2_anatase()
- structure_tio2_rutile()
"""
mutable struct Structure_Ti <: Str
    atom_vec::Matrix{Float64}
    cell_vec::Matrix{Float64}
    atom_type
    atom_name
    atom_charge
    para_mass
    num_atoms
    num_atom_types
    vec_type_id
    res_name
end

"""
    function structure_tio2_anatase()
Do this will generate a Structure_Ti type which contains all information needed to build a TiO2 moedel in structure of Anatase.
"""
function structure_tio2_anatase()
    atom_vec = [
        0.500	0.000	0.000
        0.500	0.500	0.250
        0.000	0.500	0.500
        0.000	0.000	0.750
        0.500	0.000	0.208
        0.500	0.500	0.458
        0.500	0.500	0.042
        0.500	0.000	0.792
        0.000	0.500	0.292
        0.000	0.000	0.958
        0.000	0.500	0.708
        0.000	0.000	0.542
    ]
    cell_vec = [
        3.7842      0           0
        0           3.7842      0
        0           0           9.5146
    ]
    atom_type = vcat([1 for n = 1:4], [2 for n = 5:12])
    atom_charge = 0 .* atom_type
    para_mass = [47.867, 15.999] # Ti: 47.867, O: 15.999. Unit: g/mol
    atom_name = split("Ti O")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "TIO"
    Structure_Ti(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end

"""
    function structure_tio2_rutile()
Do this will generate a Structure_Ti type which contains all information needed to build a TiO2 moedel in structure of Rutile.
"""
function structure_tio2_rutile()
    atom_vec = [
        0.804921  0.804921  0.0
        0.305029  0.305029  0.499831
        0.610059  0.610059  0.499831
        0.110168  0.499891  0.0
        0.499891  0.110168  0.0
        0.0       0.0       0.499831
    ]
    cell_vec = [
        4.593  0.0    0.0
        0.0    4.593  0.0
        0.0    0.0    2.959
    ]
    atom_type = vcat([1 for n = 1:2], [2 for n = 1:4])
    atom_charge = 0 .* atom_type
    para_mass = [47.867, 15.999] # Ti: 47.867, O: 15.999. Unit: g/mol
    atom_name = split("Ti O")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "TIO"
    Structure_Ti(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end