text:
import Base: size

using UnPack

import CartesianGrids: cellsize, origin
import RigidBodyTools: assign_velocity!
import ImmersedLayers: normals, areas

"""
    HeatConductionParameters(ρ::Float64,c::Float64,k::Float64,[thickness=nothing])

Set the density `ρ`, specific heat `c`, heat conductivity `k`, and optionally,
the thickness.
"""
struct HeatConductionParameters
    ρ :: Float64
    c :: Float64
    k :: Float64
    α :: Float64
    d :: Union{Float64,Nothing}
end
HeatConductionParameters(ρ,c,k,d) = HeatConductionParameters(ρ,c,k,k/(ρ*c),d)
HeatConductionParameters(ρ,c,k;thickness=nothing) = HeatConductionParameters(ρ,c,k,thickness)

include("operators/heaters.jl")
include("operators/linesource.jl")

"""
$(TYPEDEF)

A system type that utilizes a grid of `NX` x `NY` dual cells and `N` Lagrange forcing
points to solve the discrete transient heat conduction equation. The
parameter `static_points` specifies whether the forcing points remain static in the
grid. It should be set to `false` if a supplied motion requires that the points move.

# Constructors:

`HeatConduction(params,Δx,xlimits,ylimits,Δt
              [,pulses=nothing])` specifies the Reynolds number `Re`, the grid
              spacing `Δx`, the dimensions of the domain in the tuples `xlimits`
              and `ylimits` (excluding the ghost cells), and the time step size `Δt`.
              The other arguments are optional. The `freestream` argument can be
              passed as either a tuple (a static freestream) or a `RigidBodyMotion`
              for a time-varying freestream. The `pulses` argument can be
              used to pass in one or more spatiotemporal pulses.

`HeatConduction(Re,Δx,xlimits,ylimits,Δt,bodies::Body/BodyList
              [,problem_side=ExternalInternalProblem]
              [,ddftype=CartesianGrids.Yang3])` passes the body
              information. This constructor
              sets the motions of the body/ies to be stationary.
              The same optional arguments used for the basic constructor
              also apply for this one. In addition, the `problem_side` can be set to
              `ExternalProblem` (default), `InternalProblem`, or `ExternalInternalProblem`.
              However, it is forced to `ExternalInternalProblem` for open Bodies
              (like `Plate` type).

`HeatConduction(Re,Δx,xlimits,ylimits,Δt,bodies::Body/BodyList,
              motions::RigidBodyMotion/RigidMotionList
              [,static_points=false])`
              passes the body and associated motion information.
              The list of motions must be the same length as the list of bodies.
              The same optional arguments used for the other constructors
              also apply for this one. In addition, `static_points` can
              be set to `true` if the supplied motion should not cause the
              points to move.

"""
mutable struct HeatConduction{NX, NY, N, MT<:PointMotionType, SD<:ProblemSide, DDF<:CartesianGrids.DDFType, FT, SP} #, RKT}
    # Physical Parameters
    "Heat transfer parameters"
    params::HeatConductionParameters

    "Bodies"
    bodies::Union{BodyList,Nothing}
    "BC type"
    bctype::Type
    "Body temperatures"
    bodytemps::Union{RigidMotionList,Nothing}
    "Volumetric heat inputs"
    qforce::Union{Vector{ModulatedField},Nothing}

    "Prescribed area heat flux"
    qflux::Union{Vector{PrescribedHeatFlux},Nothing}
    "Prescribed area h*dT heat model"
    qhdT::Union{Vector{PrescribedHeatModel},Nothing}
    "Prescribed line source"
    qline::Union{Vector{PrescribedLineSource},Nothing}

    # Discretization
    "Grid metadata"
    grid::CartesianGrids.PhysicalGrid{2}
    "Time step"
    Δt::Float64

    # Layers
    dlc::Union{DoubleLayer,Nothing} # used for heat flux surface terms

    # Coordinate data, if present
    points::VectorData{N,Float64}
    normals::VectorData{N,Float64}

    # Pre-stored regularization and interpolation matrices (if present)
    Rc::Union{RegularizationMatrix,Nothing} # cell centers
    Ec::Union{InterpolationMatrix,Nothing}
    Rf::Union{RegularizationMatrix,Nothing} # cell faces
    Ef::Union{InterpolationMatrix,Nothing}
    Cc::Union{Matrix,Nothing}

    # Operators
    f :: FT

    # state vector
    state_prototype :: SP

    # Cache space
    Sc::Nodes{Primal, NX, NY,Float64}
    Sb::ScalarData{N,Float64}
    Vf::Edges{Primal, NX, NY,Float64}
    gradTf::Edges{Primal, NX, NY,Float64}
    Vb::VectorData{N,Float64}
    gradTb::VectorData{N,Float64}
    ΔTs::ScalarData{N,Float64}
    σ::ScalarData{N,Float64}
    τ::ScalarData{N,Float64}

end

function HeatConduction(params::HeatConductionParameters, Δx::Real, xlimits::Tuple{Real,Real},ylimits::Tuple{Real,Real}, Δt::Real;
                       bodies::Union{BodyList,Nothing} = nothing,
                       bodytemps::Union{RigidMotionList,Nothing} = nothing,
                       qforce::PT = nothing,
                       qflux = nothing,
                       qmodel = nothing,
                       qline = nothing,
                       static_points = true,
                       bctype = AdiabaticBC,
                       problem_side::Type{SD} = InternalProblem,
                       ddftype=CartesianGrids.Yang3) where {PT,SD<:ProblemSide}

    @unpack α = params

    g = PhysicalGrid(xlimits,ylimits,Δx)
    NX, NY = size(g)

    Fo = α*Δt/Δx^2

    Sc = Nodes{Primal,NX,NY,Float64}()
    Vf = Edges{Primal,NX,NY,Float64}()
    gradTf = Edges{Primal,NX,NY,Float64}()

    #=
    # Set up buffers
    Vf = Edges{Primal,NX,NY,Float64}()
    Vv = Edges{Primal,NX,NY,Float64}()
    Vn = Edges{Primal,NX,NY,Float64}()
    Wn = Nodes{Dual,NX,NY,Float64}()
    Vtf = EdgeGradient{Primal,Dual,NX,NY,Float64}()
    DVf = EdgeGradient{Primal,Dual,NX,NY,Float64}()
    VDVf = EdgeGradient{Primal,Dual,NX,NY,Float64}()
    =#

    #L = plan_laplacian(Sn,with_inverse=true)

    qforcefields = _process_forcing(qforce,Sc,g)

    # for now, if there are any bodies that are Open,
    # then force problem_side to ExternalInternalProblem.
    # but should be more flexible here
    problem_side_internal = _any_open_bodies(bodies) ? ExternalInternalProblem : problem_side

    N = numpts(bodies)

    Sb = ScalarData(N)
    Vb = VectorData(N)
    gradTb = VectorData(N)
    ΔTs = ScalarData(N)
    σ = ScalarData(N)
    τ = ScalarData(N)
    Cc = nothing

     points, normals, dlc, Rc, Ec, Rf, Ef, Cc =
            _immersion_operators(bodies,g,problem_side_internal,bctype,ddftype,Sc,Sb,Vf,Vb)

    conduction_L = plan_laplacian(Sc,factor=α/Δx^2)

    if isnothing(bodies)
      state_prototype = solvector(state=Sc)
      f = ConstrainedODEFunction(heatconduction_rhs!,conduction_L,_func_cache=state_prototype)
    else
      if static_points
        if bctype <: DirichletBC
          state_prototype = solvector(state=Sc,constraint=σ)
          f = ConstrainedODEFunction(heatconduction_rhs!,temperature_bc_constraint_rhs!,
                                      temperature_constraint_force!,temperature_constraint_op!,
                                      conduction_L,_func_cache=state_prototype)
        elseif bctype <: NeumannBC
          state_prototype = solvector(state=Sc)
          f = ConstrainedODEFunction(heatconduction_rhs!,conduction_L,_func_cache=state_prototype)
        end
      #=
        else
        state_prototype = solvector(state=Sc,constraint=σ,aux_state=zero_body_state(bodies))
        rhs! = ConstrainedSystems.r1vector(state_r1 = ns_rhs!,aux_r1 = rigid_body_rhs!)
        f = ConstrainedODEFunction(rhs!,temperature_bc_constraint_rhs!,
                                      temperature_constraint_force!,temperature_constraint_op!,
                                      conduction_L,_func_cache=state_prototype,
                                      param_update_func=update_immersion_operators!)
        =#
      end

    end

    HeatConduction{NX, NY, N, _motiontype(static_points), problem_side_internal, ddftype, typeof(f), typeof(state_prototype)}(
                          params, bodies, bctype, bodytemps, qforce,
                          _heat_flux_regions(qflux,Sc,g),
                          _heat_model_regions(qmodel,Sc,g),
                          _line_source(qline,Sc,g),
                          g, Δt, # rk,
                          dlc,
                          points, normals, Rc, Ec, Rf, Ef, Cc,
                          f,state_prototype,Sc,Sb,Vf,gradTf,Vb,gradTb,ΔTs,σ,τ)

end

#=
HeatConduction(Re,Δx,xlim,ylim,Δt,bodies::BodyList;
        motions=RigidMotionList(map(x -> RigidBodyMotion(0.0,0.0),bodies)),kwargs...) =
        HeatConduction(Re,Δx,xlim,ylim,Δt;bodies=bodies,motions=motions,kwargs...)

HeatConduction(Re,Δx,xlim,ylim,Δt,body::Body;kwargs...) =
        HeatConduction(Re,Δx,xlim,ylim,Δt,BodyList([body]);kwargs...)

function HeatConduction(Re,Δx,xlim,ylim,Δt,bodies::BodyList,motions::RigidMotionList;static_points=false,kwargs...)
    length(bodies) == length(motions) || error("Inconsistent lengths of bodies and motions lists")
    HeatConduction(Re,Δx,xlim,ylim,Δt,bodies;motions=motions,static_points=static_points,kwargs...)
end

HeatConduction(Re,Δx,xlim,ylim,Δt,body::Body,motion::RigidBodyMotion;static_points=false,kwargs...) =
        HeatConduction(Re,Δx,xlim,ylim,Δt,BodyList([body]),RigidMotionList([motion]);static_points=static_points,kwargs...)

=#

function Base.show(io::IO, sys::HeatConduction{NX,NY,N,MT,SD}) where {NX,NY,N,MT,SD}
    mtype = (MT == StaticPoints) ? "static" : "moving"
    sdmsg = (N == 0) ? "Unbounded" : ((SD == ExternalProblem) ? "External problem" : ((SD == InternalProblem) ? "Internal problem" : "External/internal"))
    println(io, "$sdmsg Heat conduction system on a grid of size $NX x $NY and $N $mtype immersed points")
    if N > 0
      bdmsg = (length(sys.bodies) == 1) ? "1 body" : "$(length(sys.bodies)) bodies"
      println(io, "   $bdmsg")
    end
end

# Routines to set up the immersion operators

function _immersion_operators(bodies::BodyList,g::PhysicalGrid,problem_side::Type{SD},bctype::Type{BC},ddftype,
                                Sc::ScalarGridData{NX,NY},Sb::ScalarData{N},Vf::VectorGridData{NX,NY},Vb::VectorData{N}) where {NX,NY,N,SD<:ComputationalHeatTransfer.ProblemSide,BC<:AbstractBC}

  points = VectorData(collect(bodies))
  numpts(bodies) == N || error("Inconsistent size of bodies")

  body_areas = areas(bodies)
  body_normals = normals(bodies)

  if !(problem_side==ExternalInternalProblem)
    dlc = DoubleLayer(bodies,g,Sc)
  else
    dlc = nothing
  end

  regop = _regularization(points,g,bodies,ddftype)

  Rc = RegularizationMatrix(regop,Sb,Sc) # Used by B₁ᵀ in temperature constraints
  Ec = InterpolationMatrix(regop,Sc,Sb) # Used by constraint_rhs! and B₂
  Rf = RegularizationMatrix(regop,Vb,Vf) # Used by B₁ᵀ in heat flux constraints
  Ef = InterpolationMatrix(regop,Vf,Vb) # Used by constraint_rhs! and B₂

  Cc = _neumann_matrix(bctype,body_normals,Rf,Ef,Sc,Sb,Vf,Vb)

  return points, body_normals, dlc, Rc, Ec, Rf, Ef, Cc
end

_immersion_operators(::Nothing,a...) =
    VectorData(0), VectorData(0), nothing, nothing, nothing, nothing, nothing, nothing

# For updating the system with body data

function update_immersion_operators!(sys::HeatConduction{NX,NY,N,MT,SD,DDF},bodies::BodyList) where {NX,NY,N,MT,SD,DDF<:CartesianGrids.DDFType}
    sys.bodies = deepcopy(bodies)
    sys.points, sys.dlc, sys.Rc, sys.Ec, sys.Rf, sys.Ef =
      _immersion_operators(sys.bodies,sys.grid,SD,DDF,sys.Sc,sys.Sb)
    return sys
end

update_immersion_operators!(sys::HeatConduction,body::Body) = update_immersion_operators!(sys,BodyList([body]))

function update_immersion_operators!(sys::HeatConduction,x::AbstractVector)
    tl! = RigidTransformList(x)
    tl!(sys.bodies)
    update_immersion_operators!(sys,sys.bodies)
end

# The form passed to ConstrainedODEFunction
update_immersion_operators!(sys::HeatConduction,u,sys_old::HeatConduction,t) =
    update_immersion_operators!(sys,aux_state(u))

"""
    setstepsizes(α[,gridPe=2][,cfl=0.5][,fourier=0.5]) -> Float64, Float64

Set the grid cell spacing and time step size based on the thermal diffusivity `α`,
the grid Peclet number `gridPe`, cfl number `cfl`, and grid Fourier number `fourier`.
The last three parameters all have default values.

# Example

Here is an example of setting parameters based on Reynolds number 100 (with
  default choices for grid Reynolds number, CFL number, and Fourier number):
```jldoctest
julia> Δx, Δt = setstepsizes(100)
(0.02, 0.01)
```
"""
function setstepsizes(α::Real; gridPe = 2.0, cfl = 0.5, fourier = 0.5)
    Δx = α*gridPe
    Δt = fourier*Δx^2
    return Δx, Δt
end

# some convenience functions
"""
    size(sys::HeatConduction,d::Int) -> Int

Return the number of indices of the grid used by `sys` along dimension `d`.
"""
size(sys::HeatConduction{NX,NY},d::Int) where {NX,NY} = d == 1 ? NX : NY

"""
    size(sys::HeatConduction) -> Tuple{Int,Int}

Return a tuple of the number of indices of the grid used by `sys`
"""
size(sys::HeatConduction{NX,NY}) where {NX,NY} = (size(sys,1),size(sys,2))

"""
    cellsize(sys::HeatConduction) -> Float64

Return the grid cell size of system `sys`
"""
cellsize(sys::HeatConduction) = cellsize(sys.grid)

"""
    timestep(sys::HeatConduction) -> Float64

Return the time step size of system `sys`
"""
timestep(sys::HeatConduction) = sys.Δt

"""
    origin(sys::HeatConduction) -> Tuple{Int,Int}

Return a tuple of the indices of the primal node that corresponds to the
physical origin of the coordinate system used by `sys`. Note that these
indices need not lie inside the range of indices occupied by the grid.
For example, if the range of physical coordinates occupied by the grid
is (1.0,3.0) x (2.0,4.0), then the origin is not inside the grid.
"""
origin(sys::HeatConduction) = origin(sys.grid)

"""
    timerange(tf,sys::HeatConduction)

Create a range of times, starting at the t = Δt (the time step of `sys`),
and ending at t = `tf`.
"""
timerange(tf,sys) = timestep(sys):timestep(sys):tf

# Wrap the output of the motion evaluation in VectorData
@inline assign_velocity!(u::VectorData,a...) = (assign_velocity!(u.u,u.v,a...); u)

@inline normals(sys::HeatConduction) = (!isnothing(sys.bodies)) ? normals(sys.bodies) : nothing

"""
    newstate(sys::HeatConduction)

Return a new (zero) instance of the state vector for `sys`.
"""
newstate(sys::HeatConduction) = zero(sys.state_prototype)

"""
    newstate(s::AbstractSpatialField,sys::HeatConduction)

Return an instance of the state vector for `sys`, assigned the
data in the spatial field `s`.
"""
function newstate(s::AbstractSpatialField,sys::HeatConduction)
  u = newstate(sys)
  gf = GeneratedField(state(u),s,sys.grid)
  state(u) .= cellsize(sys)*gf()
  return u
end

# Other functions
_motiontype(isstatic::Bool) = isstatic ? StaticPoints : MovingPoints

_body_closure_type(b::T) where {T<:Body{N,C}} where {N,C} = C

_any_open_bodies(nothing) = false
_any_open_bodies(bodies::BodyList) =  any(b -> _body_closure_type(b) == RigidBodyTools.OpenBody,bodies)

_regularization(sys::HeatConduction{NX, NY, N, MT, SD, DDF}) where {NX,NY,N,MT,SD,DDF} =
        _regularization(sys.points,sys.grid,sys.bodies,DDF)

_regularization(points,g,bodies,ddftype) = Regularize(points,cellsize(g),
                                I0=CartesianGrids.origin(g),weights=areas(bodies).data,ddftype=ddftype)

function _heat_flux_regions(qparams::Vector{<:PrescribedHeatFluxRegion},u::ScalarGridData,g::PhysicalGrid)
  qflux = PrescribedHeatFlux[]
  for qp in qparams
    push!(qflux,PrescribedHeatFlux(qp,u,g))
  end
  qflux
end

_heat_flux_regions(qparams::PrescribedHeatFluxRegion,u::ScalarGridData,g::PhysicalGrid) =
    _heat_flux_regions([qparams],u,g)

_heat_flux_regions(::Nothing,u,g) = nothing

function _heat_model_regions(qparams::Vector{<:PrescribedHeatModelRegion},u::ScalarGridData,g::PhysicalGrid)
  qmodel = PrescribedHeatModel[]
  for qp in qparams
    push!(qmodel,PrescribedHeatModel(qp,u,g))
  end
  qmodel
end

_heat_model_regions(qparams::PrescribedHeatModelRegion,u::ScalarGridData,g::PhysicalGrid) =
    _heat_model_regions([qparams],u,g)

_heat_model_regions(::Nothing,u,g) = nothing

function _line_source(lineparams::Vector{<:LineSourceParams},u::ScalarGridData,g::PhysicalGrid)
  qline = PrescribedLineSource[]
  for p in lineparams
    push!(qline,PrescribedLineSource(p,u,g))
  end
  qline
end

_line_source(lineparams::LineSourceParams,u::ScalarGridData,g::PhysicalGrid) =
    _line_source([lineparams],u,g)

_line_source(::Nothing,u,g) = nothing

function set_linesource_strength!(sys::HeatConduction,q::Vector{T}) where {T<:Real}
    @unpack qline = sys
    if isnothing(qline)
      return sys
    else
      total_len = mapreduce(ql -> length(ql.q),+,qline)
      @assert total_len == length(q)
      qcpy = copy(q)
      first_index = 0
      for ql in qline
        ql.q .= view(q,first_index+1:first_index+length(ql.q))
        first_index += length(ql.q)
      end
    end
    return sys
end

# Form the matrix n.Ef*Rf*n for implementing Neumann conditions
_neumann_matrix(a...) = nothing

function _neumann_matrix(::Type{<:NeumannBC},normals,Rf::RegularizationMatrix,Ef::InterpolationMatrix,
                          Sc::ScalarGridData{NX,NY},Sb::ScalarData{N},
                          Vf::VectorGridData{NX,NY},Vb::VectorData{N}) where {NX,NY,N}

  M = Matrix{Float64}(undef,N,N)

  fill!(Sb,0.0)
  for i in 1:N
    Sb[i] = 1.0
    product!(Vb,normals,Sb)
    Vf .= Rf*Vb
    Vb .= Ef*Vf
    pointwise_dot!(Sb,Vb,normals)
    M[:,i] .= Sb
    fill!(Sb,0.0)
  end
  return M
end

include("operators/surfacetemperatures.jl")
include("operators/fields.jl")
#include("operators/pointforce.jl")
include("operators/basicoperators.jl")
include("operators/surfaceoperators.jl")
include("operators/movingbodyoperators.jl")
include("operators/timemarching.jl")


text:
# Autogenerated wrapper script for Qt_jll for x86_64-linux-gnu-cxx03
export libqt53danimation, libqt53dcore, libqt53dextras, libqt53dinput, libqt53dlogic, libqt53dquick, libqt53dquickanimation, libqt53dquickextras, libqt53dquickinput, libqt53dquickrender, libqt53dquickscene2d, libqt53drender, libqt5bodymovin, libqt5charts, libqt5concurrent, libqt5core, libqt5dbus, libqt5gui, libqt5multimedia, libqt5multimediaquick, libqt5multimediawidgets, libqt5network, libqt5opengl, libqt5printsupport, libqt5qml, libqt5qmlmodels, libqt5qmlworkerscript, libqt5quick, libqt5quick3d, libqt5quick3dassetimport, libqt5quick3drender, libqt5quick3druntimerender, libqt5quick3dutils, libqt5quickcontrols2, libqt5quickparticles, libqt5quickshapes, libqt5quicktemplates2, libqt5quicktest, libqt5quickwidgets, libqt5sql, libqt5svg, libqt5test, libqt5widgets, libqt5xml, libqt5xmlpatterns

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt")
JLLWrappers.@declare_library_product(libqt53danimation, "libQt53DAnimation.so.5")
JLLWrappers.@declare_library_product(libqt53dcore, "libQt53DCore.so.5")
JLLWrappers.@declare_library_product(libqt53dextras, "libQt53DExtras.so.5")
JLLWrappers.@declare_library_product(libqt53dinput, "libQt53DInput.so.5")
JLLWrappers.@declare_library_product(libqt53dlogic, "libQt53DLogic.so.5")
JLLWrappers.@declare_library_product(libqt53dquick, "libQt53DQuick.so.5")
JLLWrappers.@declare_library_product(libqt53dquickanimation, "libQt53DQuickAnimation.so.5")
JLLWrappers.@declare_library_product(libqt53dquickextras, "libQt53DQuickExtras.so.5")
JLLWrappers.@declare_library_product(libqt53dquickinput, "libQt53DQuickInput.so.5")
JLLWrappers.@declare_library_product(libqt53dquickrender, "libQt53DQuickRender.so.5")
JLLWrappers.@declare_library_product(libqt53dquickscene2d, "libQt53DQuickScene2D.so.5")
JLLWrappers.@declare_library_product(libqt53drender, "libQt53DRender.so.5")
JLLWrappers.@declare_library_product(libqt5bodymovin, "libQt5Bodymovin.so.5")
JLLWrappers.@declare_library_product(libqt5charts, "libQt5Charts.so.5")
JLLWrappers.@declare_library_product(libqt5concurrent, "libQt5Concurrent.so.5")
JLLWrappers.@declare_library_product(libqt5core, "libQt5Core.so.5")
JLLWrappers.@declare_library_product(libqt5dbus, "libQt5DBus.so.5")
JLLWrappers.@declare_library_product(libqt5gui, "libQt5Gui.so.5")
JLLWrappers.@declare_library_product(libqt5multimedia, "libQt5Multimedia.so.5")
JLLWrappers.@declare_library_product(libqt5multimediaquick, "libQt5MultimediaQuick.so.5")
JLLWrappers.@declare_library_product(libqt5multimediawidgets, "libQt5MultimediaWidgets.so.5")
JLLWrappers.@declare_library_product(libqt5network, "libQt5Network.so.5")
JLLWrappers.@declare_library_product(libqt5opengl, "libQt5OpenGL.so.5")
JLLWrappers.@declare_library_product(libqt5printsupport, "libQt5PrintSupport.so.5")
JLLWrappers.@declare_library_product(libqt5qml, "libQt5Qml.so.5")
JLLWrappers.@declare_library_product(libqt5qmlmodels, "libQt5QmlModels.so.5")
JLLWrappers.@declare_library_product(libqt5qmlworkerscript, "libQt5QmlWorkerScript.so.5")
JLLWrappers.@declare_library_product(libqt5quick, "libQt5Quick.so.5")
JLLWrappers.@declare_library_product(libqt5quick3d, "libQt5Quick3D.so.5")
JLLWrappers.@declare_library_product(libqt5quick3dassetimport, "libQt5Quick3DAssetImport.so.5")
JLLWrappers.@declare_library_product(libqt5quick3drender, "libQt5Quick3DRender.so.5")
JLLWrappers.@declare_library_product(libqt5quick3druntimerender, "libQt5Quick3DRuntimeRender.so.5")
JLLWrappers.@declare_library_product(libqt5quick3dutils, "libQt5Quick3DUtils.so.5")
JLLWrappers.@declare_library_product(libqt5quickcontrols2, "libQt5QuickControls2.so.5")
JLLWrappers.@declare_library_product(libqt5quickparticles, "libQt5QuickParticles.so.5")
JLLWrappers.@declare_library_product(libqt5quickshapes, "libQt5QuickShapes.so.5")
JLLWrappers.@declare_library_product(libqt5quicktemplates2, "libQt5QuickTemplates2.so.5")
JLLWrappers.@declare_library_product(libqt5quicktest, "libQt5QuickTest.so.5")
JLLWrappers.@declare_library_product(libqt5quickwidgets, "libQt5QuickWidgets.so.5")
JLLWrappers.@declare_library_product(libqt5sql, "libQt5Sql.so.5")
JLLWrappers.@declare_library_product(libqt5svg, "libQt5Svg.so.5")
JLLWrappers.@declare_library_product(libqt5test, "libQt5Test.so.5")
JLLWrappers.@declare_library_product(libqt5widgets, "libQt5Widgets.so.5")
JLLWrappers.@declare_library_product(libqt5xml, "libQt5Xml.so.5")
JLLWrappers.@declare_library_product(libqt5xmlpatterns, "libQt5XmlPatterns.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt53danimation,
        "lib/libQt53DAnimation.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dcore,
        "lib/libQt53DCore.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dextras,
        "lib/libQt53DExtras.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dinput,
        "lib/libQt53DInput.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dlogic,
        "lib/libQt53DLogic.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquick,
        "lib/libQt53DQuick.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickanimation,
        "lib/libQt53DQuickAnimation.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickextras,
        "lib/libQt53DQuickExtras.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickinput,
        "lib/libQt53DQuickInput.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickrender,
        "lib/libQt53DQuickRender.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickscene2d,
        "lib/libQt53DQuickScene2D.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53drender,
        "lib/libQt53DRender.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5bodymovin,
        "lib/libQt5Bodymovin.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5charts,
        "lib/libQt5Charts.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5concurrent,
        "lib/libQt5Concurrent.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5core,
        "lib/libQt5Core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5dbus,
        "lib/libQt5DBus.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5gui,
        "lib/libQt5Gui.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimedia,
        "lib/libQt5Multimedia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimediaquick,
        "lib/libQt5MultimediaQuick.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimediawidgets,
        "lib/libQt5MultimediaWidgets.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5network,
        "lib/libQt5Network.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5opengl,
        "lib/libQt5OpenGL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5printsupport,
        "lib/libQt5PrintSupport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qml,
        "lib/libQt5Qml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlmodels,
        "lib/libQt5QmlModels.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlworkerscript,
        "lib/libQt5QmlWorkerScript.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick,
        "lib/libQt5Quick.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3d,
        "lib/libQt5Quick3D.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3dassetimport,
        "lib/libQt5Quick3DAssetImport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3drender,
        "lib/libQt5Quick3DRender.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3druntimerender,
        "lib/libQt5Quick3DRuntimeRender.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3dutils,
        "lib/libQt5Quick3DUtils.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickcontrols2,
        "lib/libQt5QuickControls2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickparticles,
        "lib/libQt5QuickParticles.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickshapes,
        "lib/libQt5QuickShapes.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktemplates2,
        "lib/libQt5QuickTemplates2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktest,
        "lib/libQt5QuickTest.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickwidgets,
        "lib/libQt5QuickWidgets.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5sql,
        "lib/libQt5Sql.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5svg,
        "lib/libQt5Svg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5test,
        "lib/libQt5Test.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5widgets,
        "lib/libQt5Widgets.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xml,
        "lib/libQt5Xml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xmlpatterns,
        "lib/libQt5XmlPatterns.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
struct FixedKernel{K<:Kernel, NFREE} <: Kernel
    kernel::K
    free::SVector{NFREE, Int}
end

@deprecate FixedKern FixedKernel

get_params(k::FixedKernel) = convert(Vector{Float64}, get_params(k.kernel)[k.free])
get_param_names(k::FixedKernel) = get_param_names(k.kernel)[k.free]
function set_params!(k::FixedKernel, hyp)
    if length(k.free) == 0
        return
    end
    p = get_params(k.kernel)
    p[k.free] = hyp
    set_params!(k.kernel, p)
end
num_params(k::FixedKernel{K,NFREE}) where {K,NFREE} = NFREE

# convenience functions to fix a parameter
function FixedKernel(k::Kernel, free::AbstractVector{<:Int})
    npars = length(free)
    sfree = SVector{npars, Int}(free)
    return FixedKernel(k, sfree)
end
function fix(k::Kernel, par::Symbol)
    npars = num_params(k)
    free = collect(1:npars)
    names = get_param_names(k)
    tofix = findfirst(==(par), names)
    tofix == nothing || deleteat!(free, tofix)
    return FixedKernel(k, free)
end

function fix(k::FixedKernel, par::Symbol)
    free_old = k.free
    names = get_param_names(k.kernel)
    free_new = [f for f in free_old if names[f] != par]
    return FixedKernel(k.kernel, free_new)
end
function fix(k::Kernel)
    return FixedKernel(k, Int[])
end

# convenience functions to free a parameter
function free(k::FixedKernel)
    return k.kernel
end
function free(k::FixedKernel, par::Symbol)
    all_names = get_param_names(k.kernel)
    ipar = findfirst(==(par), all_names)
    if ipar == nothing || ipar ∉ k.free
        free = k.free
    else
        free = sort(unique([k.free; ipar]))
    end

    return FixedKernel(k.kernel, free)
end

function grad_slice!(dK::AbstractMatrix, k::FixedKernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, p::Int)
    return grad_slice!(dK, k.kernel, X1, X2, data, k.free[p])
end
@inline function dKij_dθp(fk::FixedKernel,X1::AbstractMatrix, X2::AbstractMatrix,data::KernelData,i::Int,j::Int,p::Int,dim::Int)
    return dKij_dθp(fk.kernel, X1, X2, data, i, j, fk.free[p], dim)
end

# delegate everything else to the wrapped kernel
@inline cov_ij(fk::FixedKernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, dim::Int) = cov_ij(fk.kernel, X1, X2, data, i, j, dim)
cov(fk::FixedKernel, x::AbstractVector, y::AbstractVector) = cov(fk.kernel, x, y)
KernelData(fk::FixedKernel, X1::AbstractMatrix, X2::AbstractMatrix) = KernelData(fk.kernel, X1, X2)
kernel_data_key(fk::FixedKernel, X1::AbstractMatrix, X2::AbstractMatrix) = kernel_data_key(fk.kernel, X1, X2)

##########
# Priors #
##########

function get_priors(k::FixedKernel)
    free_priors = get_priors(k.kernel)
    if isempty(free_priors) return []
    else
        return free_priors[k.free]
    end
end

function set_priors!(k::FixedKernel, priors)
    p = get_priors(k.kernel)
    p[k.free] = priors
    set_priors!(k.kernel, p)
end

function prior_logpdf(k::FixedKernel)
    return 0.0
end

function prior_gradlogpdf(k::FixedKernel)
    return zeros(num_params(k))
end


text:
#
# Tcl.jl -
#
# Tcl interface for Julia.
#

# Mist be a bare module because it implements its own `eval` method available
# as `Tcl.eval`.
baremodule Tcl

# Only export symbols which are prefixed with `Tcl`, `TCL_`, `Tk`, `Ttk` or
# `TK_`, other "public" symbols will be available with the `Tcl.` prefix.
const __EXPORTS = (
    :TclError,
    :TclInterp,
    :TclObj,
    :TclObjList,
    :TclObjCommand,
    :TclStatus,
    :TCL_OK,
    :TCL_ERROR,
    :TCL_RETURN,
    :TCL_BREAK,
    :TCL_CONTINUE,
    :TCL_GLOBAL_ONLY,
    :TCL_NAMESPACE_ONLY,
    :TCL_APPEND_VALUE,
    :TCL_LIST_ELEMENT,
    :TCL_LEAVE_ERR_MSG,
    :TCL_DONT_WAIT,
    :TCL_WINDOW_EVENTS,
    :TCL_FILE_EVENTS,
    :TCL_TIMER_EVENTS,
    :TCL_IDLE_EVENTS,
    :TCL_ALL_EVENTS,
    :TCL_NO_EVAL,
    :TCL_EVAL_GLOBAL,
    :TCL_EVAL_DIRECT,
    :TCL_EVAL_INVOKE,
    :TCL_CANCEL_UNWIND,
    :TCL_EVAL_NOERR,
    :tkstart,
    #:@TkWidget,
    # Colors:
    :TkColor,
    :TkGray,
    :TkRGB,
    :TkBGR,
    :TkRGBA,
    :TkBGRA,
    :TkARGB,
    :TkABGR,
    # Widgets:
    :TkObject,
    :TkImage,
    :TkWidget,
    :TkRootWidget,
    :TkButton,
    :TkCanvas,
    :TkCheckbutton,
    :TkEntry,
    :TkFrame,
    :TkLabel,
    :TkLabelframe,
    :TkListbox,
    :TkMenu,
    :TkMenubutton,
    :TkMessage,
    :TkPanedwindow,
    :TkRadiobutton,
    :TkScale,
    :TkScrollbar,
    :TkSpinbox,
    :TkText,
    :TkToplevel,
    :TkWidget,
    :TtkButton,
    :TtkCheckbutton,
    :TtkCombobox,
    :TtkEntry,
    :TtkFrame,
    :TtkLabel,
    :TtkLabelframe,
    :TtkMenubutton,
    :TtkNotebook,
    :TtkPanedwindow,
    :TtkProgressbar,
    :TtkRadiobutton,
    :TtkScale,
    :TtkScrollbar,
    :TtkSeparator,
    :TtkSizegrip,
    :TtkSpinbox,
    :TtkTreeview)

# This baremodule has a special `eval` method available as `Tcl.eval`.
using Base
function eval end
function error end

if VERSION < v"0.6.0"
    # macro for raw strings (will be part of Julia 0.6, see PR #19900 at
    # https://github.com/JuliaLang/julia/pull/19900).
    export @raw_str
    macro raw_str(s); s; end
end

# Hide most implementation methods in a private module but make `Tcl` module
# available (using a relative path, in case `Tcl` is itself embedded in another
# module).
module Impl
#using ...Tcl # for public symbols
#import ...Tcl: TclInterp, TclInterpPtr, TclObj, TclObjPtr, TkImage, TclStatus
import ...Tcl

using Printf

isfile(joinpath(@__DIR__,"..","deps","deps.jl")) ||
    error("Tcl not properly installed.  Please run `Pkg.build(\"Tcl\")` to create file \"",joinpath(@__DIR__,"..","deps","deps.jl"),"\"")

include("../deps/deps.jl")
include("types.jl")
include("macros.jl")
include("calls.jl")
include("basics.jl")
include("objects.jl")
include("lists.jl")
include("variables.jl")
include("widgets.jl")
include("dialogs.jl")
include("images.jl")

end

# Import public methods, types and constants.  These will be available as
# `Tcl.$name`.
import .Impl:
    Callback,
    WideInt,
    AtomicType,
    cget,
    choosecolor,
    choosedirectory,
    colorize!,
    colorize,
    concat,
    configure,
    delete,
    deletecommand,
    exec,
    exists,
    findphoto,
    getheight,
    getinterp,
    getopenfile,
    getparent,
    getpath,
    getphotosize,
    getpixels,
    getresult,
    getsavefile,
    getvalue,
    getvar,
    getwidth,
    grid,
    isactive,
    isdeleted,
    lappend!,
    lindex,
    list,
    llength,
    messagebox,
    pack,
    place,
    resume,
    setphotosize!,
    setpixels!,
    setresult,
    setvar,
    suspend,
    threshold!,
    unsetvar

# Reexport prefixed public symbols.
for sym in __EXPORTS
    @eval begin
        import .Impl: $sym
        export $sym
    end
end

#include("shortnames.jl")

end # module


text:
using Test, WAVI

@testset "Simulation" begin
    @info "Testing simulations..."
    @testset "Simulation construction" begin
    @info "testing simulation construction"

    grid = Grid()
    bed_elevation = zeros(grid.nx, grid.ny)
    model = Model(grid = grid, bed_elevation = bed_elevation)
    simulation = Simulation(model = model, timestepping_params = TimesteppingParams(end_time = 1.))
    @test simulation isa Simulation
    end

    @testset "Timestepping" begin
    @info "testing generic timestep"

    grid = Grid()
    bed_elevation = zeros(grid.nx, grid.ny)
    model = Model(grid = grid, bed_elevation = bed_elevation)
    simulation = Simulation(model = model, timestepping_params = TimesteppingParams(end_time = 1.))
    timestep!(simulation)
    @test simulation isa Simulation
    end

end