text:
"""
    FPlane{FT} <: AbstractRotation

A parameter object for constant rotation around a vertical axis.
"""
struct FPlane{FT} <: AbstractRotation
    f :: FT
end

"""
    FPlane([FT=Float64;] f=nothing, rotation_rate=Ω_Earth, latitude=nothing)

Returns a parameter object for constant rotation at the angular frequency
`f/2`, and therefore with background vorticity `f`, around a vertical axis.
If `f` is not specified, it is calculated from `rotation_rate` and
`latitude` according to the relation `f = 2*rotation_rate*sind(latitude).

By default, `rotation_rate` is assumed to be Earth's.

Also called `FPlane`, after the "f-plane" approximation for the local effect of
a planet's rotation in a planar coordinate system tangent to the planet's surface.
"""
function FPlane(FT::DataType=Float64; f=nothing, rotation_rate=Ω_Earth, latitude=nothing)

    use_f = !isnothing(f)
    use_planet_parameters = !isnothing(latitude)

    if !xor(use_f, use_planet_parameters)
        throw(ArgumentError("Either both keywords rotation_rate and latitude must be " *
                            "specified, *or* only f must be specified."))
    end

    if use_f
        return FPlane{FT}(f)
    elseif use_planet_parameters
        return FPlane{FT}(2rotation_rate*sind(latitude))
    end
end

@inline x_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = - coriolis.f * ℑxyᶠᶜᵃ(i, j, k, grid, U.v)
@inline y_f_cross_U(i, j, k, grid, coriolis::FPlane, U) =   coriolis.f * ℑxyᶜᶠᵃ(i, j, k, grid, U.u)
@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::FPlane, U) where FT = zero(FT)

Base.show(io::IO, f_plane::FPlane{FT}) where FT =
    println(io, "FPlane{$FT}: f = ", @sprintf("%.2e", f_plane.f))

"""
    NonTraditionalFPlane{FT} <: AbstractRotation

A Coriolis implementation that facilitates non-traditional Coriolis terms in the zonal
and vertical momentum equations along with the traditional Coriolis terms.
"""
struct NonTraditionalFPlane{FT} <: AbstractRotation
    fz :: FT
    fy :: FT
end

"""
    NonTraditionalFPlane([FT=Float64;] fz=nothing, fy=nothing,
                                       rotation_rate=Ω_Earth, latitude=nothing)

Returns a parameter object for constant rotation about an axis in the `y-z` plane
with `y`- and `z`-components `fy/2` and `fz/2`, and the background vorticity is
`(0, fy, fz)`.

In oceanography `fz` and `fy` represent the components of planetary voriticity which
are perpendicular and parallel to the ocean surface in a domain in which `x, y, z`
correspond to the directions east, north, and up.

If `fz` and `fy` are not specified, they are calculated from `rotation_rate` and `latitude`
according to the relations `fz = 2*rotation_rate*sind(latitude)` and
`fy = 2*rotation_rate*cosd(latitude)`, respectively. By default, `rotation_rate`
is assumed to be Earth's.
"""
function NonTraditionalFPlane(FT=Float64; fz=nothing, fy=nothing, rotation_rate=Ω_Earth, latitude=nothing)

    use_f = !isnothing(fz) && !isnothing(fy) && isnothing(latitude)
    use_planet_parameters = !isnothing(latitude) && isnothing(fz) && isnothing(fy)

    if (use_planet_parameters && use_f) || (!use_planet_parameters && !use_f)
        throw(ArgumentError("Either both rotation_rate and latitude must be " *
                            "specified, *or* both f and f′ must be specified."))
    elseif use_planet_parameters
        fz = 2rotation_rate*sind(latitude)
        fy = 2rotation_rate*cosd(latitude)
    end
    return NonTraditionalFPlane{FT}(fz, fy)
end

@inline fᶻv_minus_fʸw(i, j, k, grid, coriolis::NonTraditionalFPlane, U) = coriolis.fy * ℑzᵃᵃᶜ(i, j, k, grid, U.w) - coriolis.fz * ℑyᵃᶜᵃ(i, j, k, grid, U.v)

@inline x_f_cross_U(i, j, k, grid, coriolis::NonTraditionalFPlane, U) =   ℑxᶠᵃᵃ(i, j, k, grid, fᶻv_minus_fʸw, coriolis, U)
@inline y_f_cross_U(i, j, k, grid, coriolis::NonTraditionalFPlane, U) =   coriolis.fz * ℑxyᶜᶠᵃ(i, j, k, grid, U.u)
@inline z_f_cross_U(i, j, k, grid, coriolis::NonTraditionalFPlane, U) = - coriolis.fy * ℑxzᶜᵃᶠ(i, j, k, grid, U.u)

Base.show(io::IO, f_plane::NonTraditionalFPlane{FT}) where FT =
    println(io, "Non-Traditional FPlane{$FT}:", @sprintf("fz = %.2e, fy = %.2e", f_plane.fz, f_plane.fy))


text:
#!/usr/bin/env julia

using Documenter
using DataKnots

# Setup for doctests embedded in docstrings.
DocMeta.setdocmeta!(DataKnots, :DocTestSetup, :(using DataKnots))

# Highlight indented code blocks as Julia code.
using Documenter.Expanders: ExpanderPipeline, Selectors, Markdown, iscode
abstract type DefaultLanguage <: ExpanderPipeline end
Selectors.order(::Type{DefaultLanguage}) = 99.0
Selectors.matcher(::Type{DefaultLanguage}, node, page, doc) =
    iscode(node, "")
Selectors.runner(::Type{DefaultLanguage}, node, page, doc) =
    page.mapping[node] = Markdown.Code("julia", node.code)

makedocs(
    sitename = "DataKnots.jl",
    format = Documenter.HTML(prettyurls=(get(ENV, "CI", nothing) == "true")),
    pages = [
        "Home" => "index.md",
        "Queries for Data Analysts" => "overview.md",
        "Tutorials" => [
            "highlypaid.md"
        ],
        "Reference Manual" => [
            "primer.md",
            "tutorial.md",
            "reference.md",
        ],
        "Implementer's Guide" => [
            "vectors.md",
            "pipelines.md",
            "shapes.md",
            "knots.md",
            "queries.md",
        ],
    ],
    modules = [DataKnots])

deploydocs(
    repo = "github.com/MechanicalRabbit/DataKnots.jl.git",
)

ext: jl


text:
using PyFrom
using Test

@testset "PyFrom" begin

    # Correct syntax
    @pyfrom collections.abc import Callable
    @test Callable !== nothing

    # Correct syntax
    @pyfrom math import inf as py∞, pi as pyπ, tau
    @test pyπ ≈ π
    @test py∞ ≈ Inf
    @test tau ≈ 2π

    # Incorrect syntax
    @test_throws LoadError @macroexpand @pyfrom math import a.b.c as d
    @test_throws PyFrom.ImportPhraseException PyFrom.pyfrom(:math, :(import a.b.c as d))

    # Correct syntax, but wrong imports
    @test PyFrom.pyfrom(:math, :(import a as aa, b as bb, c)) isa Expr
    @test @macroexpand( @pyfrom math import a as aa, b as bb, c ) isa Expr
    @test_throws Exception @pyfrom math import a as aa, b as bb, c

    # Do we understand a.b.c.d correctly?
    @test PyFrom.dotpath_expansion(:(a.b.c.d)) isa Tuple
    @test PyFrom.dotpath_expansion(:(a)) isa Tuple
    @test_throws PyFrom.DotExpansionException PyFrom.dotpath_expansion(:(a."b".c.d))
    @test_throws PyFrom.DotExpansionException PyFrom.dotpath_expansion(:(a.c.d + 5))

    # Do we understand import ... correctly?
    @test PyFrom.importphrase_to_mapping(:(import a, b, c)) == [:a=>:a, :b=>:b, :c=>:c]
    @test PyFrom.importphrase_to_mapping(:(import a as A, b, c)) == [:a=>:A, :b=>:b, :c=>:c]
    @test_throws PyFrom.ImportPhraseException PyFrom.importphrase_to_mapping(:(import a as A, b.c, c))

end;


text:
"""
    SeisOMP(in,operator,param,iter,lp,imax)

Robust Orthognal Matching Pursuit algorithm with the coefficient updated by the conjugate gradient pursuit

# Arguments
- `dn`: Noisy data
- `operator`: Any operator passed the dot test
- `param`: parameters for the operator
- `iter`: Max iterations for the matching pursuit
- `lp`: lp-norm space inner product
- `imax`: Max interation for the lp-norm inner product

"""
function SeisOMP(dn,operator,param,iter,lp,imax)
    m=operator(dn,param,-1,lp,imax);
    r=dn;
    n=1;
    ind_r=[];
    ind_c=[];
    v=[];
    d=[];
    dr=zeros(size(dn))
    for i=1:iter
        g=zeros(size(m))
        gn=operator(r,param,-1,lp,imax)
        ind_ri=Tuple(findmax(abs.(gn))[2])[1]
        ind_ci=Tuple(findmax(abs.(gn))[2])[2]
        push!(ind_r,ind_ri)
        push!(ind_c,ind_ci)
        for j=length(ind_r)
            g[ind_r[j],ind_c[j]]=gn[ind_r[j],ind_c[j]]
        end
        if n==1
            d=g;
            v=operator(d,param,1,lp,imax);
        else
            w=operator(g,param,1,lp,imax)
            u=-LP_norm(vec(v),vec(w),lp,imax,1)
            d=g+u*d
            v=w+u*v
        end
        beta=LP_norm(vec(v),vec(r),lp,imax,1)
        dr=dr+beta*v
        r=r-beta*v
        n=n+1
    end
    return dr
end

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
import Base: convert

convert(::Type{ArrayExpr}, val::Number) = LazyArray(val)
convert(::Type{ArrayExpr}, val::AbstractArray) = LazyArray(val)
convert(::Type{ArrayExpr}, val::ArrayExpr) = val

"""
Converts a Julia array into an object that can be used in array expressions.
This happens automatically when an array is used in an expression with an `ArrayExpr` object,
but it can be useful to wrap it explicitly to delay the computation.
"""
lazy_array(val::Number) = convert(ArrayExpr, val)
lazy_array(val::AbstractArray) = convert(ArrayExpr, val)

function unknown_array(; shape=nothing, eltype=nothing)
    Unknown(ArrayType(shape, eltype))
end

unknown_random_state() = Unknown(RandomStateType())


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LaMEM"
version = v"1.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://bitbucket.org/bkaus/lamem.git",
    "734724d6686c203f6e0012c1f087181de3da4f53")
]

# Bash recipe for building across all platforms
script = raw"""

# Create required directories
mkdir $WORKSPACE/srcdir/lamem/bin
mkdir $WORKSPACE/srcdir/lamem/bin/opt
mkdir $WORKSPACE/srcdir/lamem/dep
mkdir $WORKSPACE/srcdir/lamem/dep/opt
mkdir $WORKSPACE/srcdir/lamem/lib
mkdir $WORKSPACE/srcdir/lamem/lib/opt

cd $WORKSPACE/srcdir/lamem/src
export PETSC_OPT=${libdir}/petsc/double_real_Int32/
make mode=opt all -j${nproc}

# compile dynamic library
make mode=opt dylib -j${nproc}

cd  $WORKSPACE/srcdir/lamem/bin/opt

# On some windows versions it automatically puts the .exe extension; on others not.
# This deals with that
if [[ -f LaMEM ]]
then
    mv LaMEM LaMEM${exeext}
fi

cp LaMEM${exeext} $WORKSPACE/srcdir/lamem/
cp LaMEM${exeext} $WORKSPACE/srcdir
cd $WORKSPACE/srcdir/lamem

# Install binaries
install -Dvm 755 LaMEM* "${bindir}/LaMEM${exeext}"
install -vm 644 src/*.h "${includedir}"
install -Dvm 755 lib/opt/LaMEMLib.dylib "${libdir}/LaMEMLib.${dlext}"

# Install license
install_license LICENSE

exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms(exclude=[Platform("i686", "windows"),
                                                                  Platform("i686", "linux"; libc = "musl")]))

# The products that we will ensure are always built
products = [
    ExecutableProduct("LaMEM", :LaMEM)
    LibraryProduct("LaMEMLib", :LaMEMLib)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="PETSc_jll", uuid="8fa3689e-f0b9-5420-9873-adf6ccf46f2d")),
    Dependency("MPICH_jll"; platforms=filter(!Sys.iswindows, platforms)),
    Dependency("MicrosoftMPI_jll"; platforms=filter(Sys.iswindows, platforms)),
    Dependency("CompilerSupportLibraries_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"10.2.0")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
#using Statistics

function stepOcean_MLDynamics!(
    ocn   :: Ocean;
    cfgs...
)

    # Unpacking
    do_qflx       = cfgs[:do_qflx]
    use_h_ML      = cfgs[:use_h_ML]
    Δt            = cfgs[:Δt]
    do_convadjust = cfgs[:do_convadjust]
    rad_scheme    = cfgs[:rad_scheme]

    ifrac   = ocn.in_flds.ifrac

    taux    = ocn.in_flds.taux
    tauy    = ocn.in_flds.tauy
    fric_u  = ocn.fric_u

    swflx   = ocn.in_flds.swflx
    nswflx  = ocn.in_flds.nswflx
    frwflx  = ocn.in_flds.frwflx
    vsflx   = ocn.in_flds.vsflx
    qflx_T  = ocn.in_flds.qflx_T
    qflx_S  = ocn.in_flds.qflx_S

#=
    println("get in hz ##### Ts: ", ocn.Ts[1:5, 48, 89])
    println("get in hz ##### Ss: ", ocn.Ss[1:5, 48, 89])
    println("get in hz ##### bs: ", ocn.bs[1:5, 48, 89])
=#

    # It is assumed here that buoyancy has already been updated.
    @loop_hor ocn i j let

        zs = ocn.cols.zs[i, j]
        Nz = ocn.Nz[i, j]

        fric_u[i, j] = √( √(taux[i, j]^2.0 + tauy[i, j]^2.0) / HOOM.ρ_sw)
        weighted_fric_u = fric_u[i, j] #* (1.0 - ifrac[i, j])

        # Pseudo code
        # Current using only Euler forward scheme:
        # 1. Determine h at t+Δt
        # 2. Determine how many layers are going to be
        #    taken away by ML.
        # 3. Cal b at t+Δt for both ML and DO
        # 4. Detect if it is buoyantly stable.
        #    Correct it (i.e. convection) if it is not.
        # 5. If convection happens, redetermine h.

        # p.s.: Need to examine carefully about the
        #       conservation of buoyancy in water column

        old_FLDO = ocn.FLDO[i, j]
        old_h_ML = ocn.h_ML[i, j]
        old_T_ML = ocn.T_ML[i, j]
        old_S_ML = ocn.S_ML[i, j]

        α = TS2α(old_T_ML, old_S_ML)
        β = TS2β(old_T_ML, old_S_ML)

        surf_Tnswflx = nswflx[i, j] / ρc_sw
        surf_Tswflx  = swflx[i, j] / ρc_sw
        surf_Jflx    = g * α * surf_Tswflx
        #surf_Sflx    = - frwflx[i, j] * ocn.S_ML[i, j] / ρ_fw
        surf_Sflx    = vsflx[i, j]
        surf_bflx    = g * ( α * surf_Tnswflx - β * surf_Sflx )

        ocn.SFLUX_top[i, j] = surf_Sflx

        new_h_ML = old_h_ML

        if use_h_ML # h_ML is datastream

            new_h_ML = ocn.in_flds.h_ML[i, j]

        else        # h_ML is prognostic

#            target_z = max( - old_h_ML - 30.0,  - ocn.h_ML_max[i, j])
#            avg_D = - old_h_ML - target_z

#=
            if (i, j) == (48, 89)
                println("before delta b ##### Ts: ", ocn.Ts[1:5, i, j])
                println("before delta b ##### Ss: ", ocn.Ss[1:5, i, j])
                println("before delta b ##### bs: ", ocn.bs[1:5, i, j])
            end
=#

#            Δb = ( (avg_D > 0.0) ? ocn.b_ML[i, j] - (
#                  OC_getIntegratedBuoyancy(ocn, i, j; target_z =   target_z)
#                - OC_getIntegratedBuoyancy(ocn, i, j; target_z = - old_h_ML)
#            ) / avg_D
#            : 0.0 )

            Δb = (old_FLDO == -1 ) ? 0.0 : ocn.b_ML[i, j] - ocn.bs[old_FLDO, i, j]

            # After convective adjustment, there still might
            # be some numerical error making Δb slightly negative
            # (the one I got is like -1e-15 ~ -1e-8). So I set a
            # tolarence δb = 3e-6 ( 0.001 K => 3e-6 m/s^2 ).
            if -3e-6 < Δb < 0.0
                Δb = 0.0
            end

            #=
            if Δb < -3e-6
                println(format("[MLDynamics] At {:d}, {:d}: Δb = {:f}. FLDO = {:d}", i, j, Δb, old_FLDO))
                println("See if T_ML and Ts are inconsistent: ")
                println("T_ML: ", ocn.T_ML[i, j])
                println("Ts  : ", ocn.Ts[:, i, j])
                println(format("ΔT  : {:e}", ocn.T_ML[i, j] - ocn.Ts[1, i, j]))
                println("See if S_ML and Ss are inconsistent: ")
                println("S_ML: ", ocn.S_ML[i, j])
                println("Ss  : ", ocn.Ss[:, i, j])
                println(format("ΔS  : {:e}", ocn.S_ML[i, j] - ocn.Ss[1, i, j]))
            end
=#
#            if Δb < 0.0
#                FLDO = ocn.FLDO[i, j]

#=
                println(format("({:d},{:d}) Averge sense Δb={:f}", i, j, Δb))
                println(format("({:d},{:d}) Jump sense Δb={:f}", i, j, (FLDO != -1) ? ocn.b_ML[i, j] - ocn.bs[FLDO, i, j] : 999 ))
                println("##### Ts: ", ocn.Ts[:, i, j])
                println("##### Ss: ", ocn.Ss[:, i, j])
                println("##### bs: ", ocn.bs[:, i, j])
=#
#            end

            new_h_ML, ocn.h_MO[i, j] = calNewMLD(;
                h_ML   = old_h_ML,
                Bf     = surf_bflx + surf_Jflx * ocn.R,
                J0     = surf_Jflx * (1.0 - ocn.R),
                fric_u = weighted_fric_u,
                Δb     = Δb,
                f      = ocn.fs[i, j],
                Δt     = Δt,
                ζ      = ocn.ζ,
                h_max  = ocn.h_ML_max[i, j],
                we_max = ocn.we_max,
            )

        end

        new_h_ML = boundMLD(new_h_ML; h_ML_max=ocn.h_ML_max[i, j], h_ML_min=ocn.h_ML_min[i, j])

        # ML
        #      i: Calculate integrated buoyancy that should
        #         be conserved purely through entrainment
        #     ii: Add to total buoyancy

        # If new_h_ML < old_h_ML, then the FLDO layer should get extra T or S due to mixing

        if new_h_ML < old_h_ML

            new_FLDO = getFLDO(zs=zs, h_ML=new_h_ML, Nz=Nz)

            if old_FLDO == -1

                # Mixing does not happen because FLDO does not exist in this case
                ocn.Ts[new_FLDO:Nz, i, j] .= ocn.T_ML[i, j]
                ocn.Ss[new_FLDO:Nz, i, j] .= ocn.S_ML[i, j]

            else
                FLDO_Δz =  -old_h_ML - zs[old_FLDO+1]
                retreat_Δz =  old_h_ML - ( (new_FLDO == old_FLDO) ? new_h_ML : (-zs[old_FLDO]) )

                ocn.Ts[old_FLDO, i, j] = (
                    ocn.Ts[old_FLDO, i, j] * FLDO_Δz + ocn.T_ML[i, j] * retreat_Δz
                ) / (FLDO_Δz + retreat_Δz)

                ocn.Ss[old_FLDO, i, j] = (
                    ocn.Ss[old_FLDO, i, j] * FLDO_Δz + ocn.S_ML[i, j] * retreat_Δz
                ) / (FLDO_Δz + retreat_Δz)
            end
        end

        if_entrainment = new_h_ML > old_h_ML

        # Calculate the effect of entrainment on SSS
        new_int_S_ML = OC_getIntegratedSalinity(   ocn, i, j; target_z = -new_h_ML)
        new_S_ML = new_int_S_ML / new_h_ML
        ocn.dSdt_ent[i, j] = (if_entrainment) ? (new_S_ML - old_S_ML) / Δt : 0.0

        # Add in external surface flux effect on SSS
        new_S_ML = (new_int_S_ML - surf_Sflx * Δt) / new_h_ML

        # Calculate the effect of entrainment on SST
        new_int_T_ML = OC_getIntegratedTemperature(ocn, i, j; target_z = -new_h_ML)
        new_T_ML = new_int_T_ML / new_h_ML
        ocn.dTdt_ent[i, j] = (if_entrainment) ? (new_T_ML - old_T_ML) / Δt : 0.0

        # Add in external surface flux effect on SST. Shortwave radiation is not included yet
        new_T_ML = (new_int_T_ML - surf_Tnswflx * Δt) / new_h_ML

        # Q-flux
        if do_qflx

            new_T_ML += qflx_T[i, j] * Δt / (ρc_sw * new_h_ML)
            new_S_ML += qflx_S[i, j] * Δt / new_h_ML

        end

        # Update mixed-layer
        OC_setMixedLayer!(
            ocn, i, j;
            T_ML=new_T_ML,
            S_ML=new_S_ML,
            h_ML=new_h_ML,
        )

#            if ocn.FLDO[i, j] > 1 && ocn.T_ML[i, j] != ocn.Ts[1, i, j]
#                println(format("UPDATE ML ERROR: ({},{}) has T_ML={:f} but Ts[1]={:f}", i, j, ocn.T_ML[i,j], ocn.Ts[1,i,j]))
#            end

        # Shortwave radiation
        if rad_scheme == :exponential
            FLDO = ocn.FLDO[i, j]
            ocn.T_ML[i, j] += - ocn.R * surf_Tswflx * Δt / new_h_ML
            ocn.Ts[1:((FLDO == -1) ? Nz : FLDO-1 ), i, j] .= ocn.T_ML[i, j]
            OC_doShortwaveRadiation!(ocn, i, j; Tswflx=(1.0 - ocn.R) * surf_Tswflx, Δt=Δt)

#            if ocn.FLDO[i, j] > 1 && ocn.T_ML[i, j] != ocn.Ts[1, i, j]
#                println(format("RADIATION ERROR: ({},{}) has T_ML={:f} but Ts[1]={:f}", i, j, ocn.T_ML[i,j], ocn.Ts[1,i,j]))
#            end
        elseif rad_scheme == :step
            FLDO = ocn.FLDO[i, j]
            ocn.T_ML[i, j] += - surf_Tswflx * Δt / new_h_ML
            ocn.Ts[1:((FLDO == -1) ? Nz : FLDO-1 ), i, j] .= ocn.T_ML[i, j]
        end

        OC_updateB!(ocn, i, j)

        if do_convadjust
            OC_doConvectiveAdjustment!(ocn, i, j;)

#            if FLDO > 1 && ocn.T_ML[i, j] != ocn.Ts[1, i, j]
#                println(format("CONV ERROR: ({},{}) has T_ML={:f} but Ts[1]={:f}", i, j, ocn.T_ML[i,j], ocn.Ts[1,i,j]))
#            end

        end

    end

end

function stepOcean_slowprocesses!(
    ocn :: Ocean;
    cfgs...
)

    Δt            = cfgs[:Δt]
    do_vert_diff  = cfgs[:do_vert_diff]
    do_horz_diff  = cfgs[:do_horz_diff]
    do_relaxation = cfgs[:do_relaxation]
    do_convadjust = cfgs[:do_convadjust]

    # Climatology relaxation
    if do_relaxation
        start_layer = ocn.deep_ocn_correction_start_layer
        @loop_hor ocn i j let
            ocn.TSAS_clim[i, j] = OC_doDeepOcnCorrectionOfT!(ocn, i, j; Δt=Δt, start_layer=start_layer)
            ocn.SSAS_clim[i, j] = OC_doDeepOcnCorrectionOfS!(ocn, i, j; Δt=Δt, start_layer=start_layer)
        end
    end

    # Vertical diffusion
    if do_vert_diff
        @loop_hor ocn i j let
            OC_doDiffusion_EulerBackward!(ocn, i, j; Δt=Δt)
        end
    end

    if do_relaxation || do_vert_diff || do_horz_diff

        @loop_hor ocn i j let
            OC_updateB!(ocn, i, j)
        end

        if do_convadjust
            @loop_hor ocn i j let
                OC_doConvectiveAdjustment!(ocn, i, j;)
            end
        end

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function axdoubling!(ax::Axis)
    #WRONG: sort!(append!(ax.ticks, ax.ticks[1:end - 1] + diff(ax.ticks) / 2); alg = QuickSort)
    #the axis could be in a non-sorted form

    doubleindex=(2:2*length(ax.ticks)).÷2
    doubleindex[2:2:end] .+=length(ax.ticks)
    append!(ax.ticks, ax.ticks[1:end - 1] + diff(ax.ticks) / 2)
    ax.ticks[1:end].=ax.ticks[doubleindex];
end

"""
    axesextend!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, axisnumber::Integer; prepend::AbstractArray=[], append::AbstractArray=[])

Extend parameter space of the selected Axis by prepend and append the grid.

Note, that the resolution and the monotonically is not checked.

# Examples
extending the second parameter of the mdbm problem:
```jldoctest
julia> axesextend!(mymdbm,2,prepend=-6.2:0.05:-2.2,append=2.2:0.2:3);
```
"""
function axesextend!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, axisnumber::Integer; prepend::AbstractArray = [], append::AbstractArray = []) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    preplength = length(prepend);
    if preplength > 0
        prepend!(mdbm.axes[axisnumber].ticks, prepend)
        for nc in mdbm.ncubes
            nc.corner[axisnumber] += preplength;
        end
    end
    if length(append) > 0
        append!(mdbm.axes[axisnumber].ticks, append)
    end
end

# function corner(nc::NCube{IT,FT,N}, T01)::Vector{MVector} where IT where FT where N
function corner(nc::NCube{IT,FT,N}, T01) where IT where FT where N
    [nc.corner .+ nc.size .* T for T in T01]
end

# function corner(ncubes::Vector{NCube{IT,FT,N}}, T01)::Vector{Vector{MVector}} where IT where FT where N
function corner(ncubes::Vector{NCube{IT,FT,N}}, T01) where IT where FT where N
    [corner(nc, T01) for nc in ncubes]
    # [nc.corner .+ nc.size .* T for nc in ncubes for T in T01]
end

# function corner(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})::Vector{Vector{SVector}}  where fcT where N where Nf where Nc
function corner(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    [corner(nc, mdbm.T01) for nc in mdbm.ncubes]
end

function getcornerval(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}#get it for all
    #getcornerval.(mdbm.ncubes,Ref(mdbm))
    map(x->((mdbm.fc) ∘ (mdbm.axes))(x...), Base.Iterators.flatten(corner(mdbm)))
end

function getcornerval(ncubes::NCube{IT,FT,N}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    #PostionsVectors=(mdbm.axes.(corner(nc,mdbm.T01)))
    # (mdbm.fc).( (mdbm.axes).(corner(nc,mdbm.T01)) )
    map(x->((mdbm.fc) ∘ (mdbm.axes))(x...), corner(ncubes, mdbm.T01))
end

function T01maker(valk::Val{kdim}) where {kdim}
    SVector{2^kdim}([
        SVector{kdim}([isodd(x ÷ (2^y)) for y in 0:(kdim - 1)])
         for  x in 0:(2^kdim - 1)])
end

function T11pinvmaker(valk::Val{kdim}) where {kdim}
    T11pinv = ([isodd(x ÷ (2^y)) for y in 0:kdim , x in 0:(2^kdim - 1)] * 2.0 .- 1.0) / (2^kdim)
    SMatrix{kdim + 1,2^kdim}(T11pinv)
end

function issingchange(FunTupleVector::AbstractVector, Nf::Integer, Nc::Integer)::Bool
    all([
        [
        any((c)->!isless(c[1][fi], zero(c[1][fi])), FunTupleVector)
        for fi in 1:Nf#length(FunTupleVector[1][1])
        ]#any positive (or zeros) value (!isless)

        [
        any((c)->!isless(zero(c[1][fi]), c[1][fi]), FunTupleVector)
        for fi in 1:Nf#length(FunTupleVector[1][1])
        ]#anany negative (or zero) value

        [
        any((c)->!isless(c[2][fi], zero(c[2][fi])), FunTupleVector)
        for fi in 1:Nc#length(FunTupleVector[1][2])
        ]#any positive (or zeros) value (!isless)
    ])#check for all condition
end

function _interpolate!(ncubes::Vector{NCube{IT,FT,N}}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, ::Type{Val{0}})  where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    isbracketing = map(nc->issingchange(getcornerval(nc, mdbm), Nf, Nc), ncubes) #TODO: parallelize
    deleteat!(ncubes,.!isbracketing)#removeing the non-bracketing ncubes

    for nc in ncubes
        nc.posinterp[:] .= zero(FT)
    end
    return nothing
end

function _interpolate!(ncubes::Vector{NCube{IT,FT,N}}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, ::Type{Val{1}})  where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    for nc in ncubes
        FunTupleVector = getcornerval(nc, mdbm)
        if Nc==0 || all(
            any((c)->!isless(c[2][fi], zero(c[2][fi])), FunTupleVector)
            for fi in 1:length(FunTupleVector[1][2])
            )# check the constraint: do wh have to compute at all?!?

            As = zero(MVector{Nf + Nc,FT})
            ns = zero(MMatrix{N,Nf + Nc,FT})

            #for f---------------------
            for kf = 1:Nf#length(FunTupleVector[1][1])
                solloc = mdbm.T11pinv * [FunTupleVector[kcubecorner][1][kf] for kcubecorner = 1:length(FunTupleVector)]
                As[kf] = solloc[end];
                ns[:,kf] .= solloc[1:end - 1];
            end

            #for c---if needed: that is, it is close to the boundary--------
            activeCostraint = 0;
            for kf = 1:Nc#length(FunTupleVector[1][2])
                #TODO: mivan a többszörös C teljesülése esetén!?!??!# ISSUE
                if any((c)->!isless(zero(c[2][kf]), c[2][kf]), FunTupleVector) && length(As) < N  # use a constraint till it reduces the dimension to zero (point) and no further
                    solloc = mdbm.T11pinv * [FunTupleVector[kcubecorner][2][kf] for kcubecorner = 1:length(FunTupleVector)]
                    activeCostraint += 1;
                    As[Nf + activeCostraint] = solloc[end];
                    ns[:,Nf + activeCostraint] .= solloc[1:end - 1];
                end
            end

            if (Nf + activeCostraint) == 0
                nc.posinterp[:] .= zero(FT)
            else
                # first two leads to error in the fucntion of constant due to the behaviour of pinv (eg. [1 0;0 1;0 0]/[1 0; 0 0; 0 0])
                # nc.posinterp[:] .= transpose(view(ns, :, 1:(Nf + activeCostraint))) \ view(As, 1:(Nf + activeCostraint));#TEST
                # nc.posinterp[:] .= transpose(ns[:, 1:(Nf + activeCostraint)]) \ As[1:(Nf + activeCostraint)];#TEST
                # nc.posinterp[:] .= ns[:, 1:(Nf + activeCostraint)] * ((transpose(ns[:, 1:(Nf + activeCostraint)]) * ns[:, 1:(Nf + activeCostraint)]) \ view(As, 1:(Nf + activeCostraint)));

                a = ns[:, 1:(Nf + activeCostraint)]
                d = view(As, 1:(Nf + activeCostraint))
                A = transpose(a) * a
                if rank(A) == minimum(size(A))
                    nc.posinterp[:] .= a * (A \ d)
                else
                    nc.posinterp[:] .= 1000.0
                end
            end
            #for c---------------------
        else
            nc.posinterp[:] .= 1000.0;#put it outside the cube!
        end
    end

    #TODO: what if it falls outside of the n-cube, it should be removed ->what shall I do with the bracketing cubes?
    # Let the user define it
     # filter!((nc)->sum((abs.(nc.posinterp)).^10.0)<(1.5 ^10.0),mdbm.ncubes)#1e6~=(2.0 ^20.0)
    normp = 5.0
    ncubetolerance = 0.7
    filter!((nc)->norm(nc.posinterp, normp) < 1.0 + ncubetolerance, mdbm.ncubes)
     #filter!((nc)->!any(isnan.(nc.posinterp)),mdbm.ncubes)

    return nothing
end

function _interpolate!(ncubes::Vector{NCube{IT,FT,N}}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, ::Type{Val{Ninterp}}) where {fcT, IT, FT, N, Ninterp, Nf, Nc, t01T, t11T, aT}
    error("order $(Ninterp) interpolation is not supperted (yet)")
end

"""
    interpolate!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}; interpolationorder::Int=1)

Interpolate the solution within the n-cube with order `interpolationorder`.
 - `interpolationorder=0` provide the midpoint of the n-cube
 - `interpolationorder=1` preform a linear fit and provide closest solution point to the mindpoint of the n-cube
"""
function interpolate!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}; interpolationorder::Int = 1) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    _interpolate!(mdbm.ncubes, mdbm, Val{interpolationorder})
end

"""
    refine!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}; directions::Vector{T}=collect(Int64,1:N)) where N where Nf where Nc where T<:Integer

Double the resolution of the axes along the provided `directions` then halving the `ncubes` accordingly.

# Examples
```jldoctest
julia> refine!(mymdbm)
julia> refine!(mymdbm,[1,1,2])
```
"""
function refine!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}; directions::Vector{T} = collect(Int64, 1:N))  where {T<:Integer, fcT, IT, FT, N, Ninterp, Nf, Nc, t01T, t11T, aT}
    doubling!(mdbm, directions)
    refinencubes!(mdbm.ncubes, directions)
    return nothing
end

"""
    doubling!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, directions::Vector{T}) where T<:Integer  where N where Nf where Nc

Double the resolution of the axes along the provided `directions`, then set the new size of the n-cubes.

# Examples
```jldoctest
julia> doubling!(mymdbm)
julia> doubling!(mymdbm,[1,2])
```
"""
function doubling!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, directions::Vector{T}) where T <: Integer where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    axdoubling!.(mdbm.axes[directions])
    for nc in mdbm.ncubes
        for dir in directions
            nc.corner[dir] = (nc.corner[dir] - 1) * 2 + 1
            nc.size[dir] *= 2
        end
    end
end

"""
     refinencubes!(ncubes::Vector{NCube{IT,FT,N}}, directions::Vector{T})

Halving the `ncubes` along the provided `directions`.

# Examples
```jldoctest
julia> refinencubes!(mymdbm.ncubes)
julia> refinencubes!(mymdbm.ncubes,mymdbm,[1,2])
```
"""
function refinencubes!(ncubes::Vector{NCube{IT,FT,N}}, directions::Vector{T}) where IT where FT where N where T <: Integer where t01T #{IT,FT,N} where IT where FT where N
    for dir in directions
        for nc in ncubes
            nc.size[dir] /= 2
        end
        NumofNCubes = length(ncubes)
        append!(ncubes, deepcopy(ncubes))
        for nc in ncubes[1:NumofNCubes]
            nc.corner[dir] = nc.corner[dir] + nc.size[dir]
        end
    end
    sort!(ncubes; alg = QuickSort)
    return nothing
end

"""
    getinterpolatedsolution(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})

Provide the interpolated coordinates of the all the detected solution (approximately where foo(x,y) == 0 and c(x,y)>0).
"""
function getinterpolatedsolution(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
[
    [
        (typeof(mdbm.axes[i].ticks).parameters[1])((mdbm.axes[i].ticks[nc.corner[i]] * (1.0 - (nc.posinterp[i] + 1.0) / 2.0) +
        mdbm.axes[i].ticks[nc.corner[i] + 1] * ((nc.posinterp[i] + 1.0) / 2.0)))
    for nc in mdbm.ncubes]#::Vector{typeof(mdbm.axes[i].ticks).parameters[1]}
for i in 1:length(mdbm.axes)]
end

"""
    getinterpolatedsolution(ncubes::Vector{NCube{IT,FT,N}},mdbm::MDBM_Problem{N,Nf,Nc})

Provide the interpolated coordinates of the detected solution for the selected n-cubes (approximately where foo(x,y) == 0 and c(x,y)>0).
"""
function getinterpolatedsolution(ncubes::Vector{NCube{IT,FT,N}}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
[
    [
        (typeof(mdbm.axes[i].ticks).parameters[1])(
        (mdbm.axes[i].ticks[nc.corner[i]]*(1.0-(nc.posinterp[i]+1.0)/2.0)+
        mdbm.axes[i].ticks[nc.corner[i]+1]*((nc.posinterp[i]+1.0)/2.0))
        )
    for nc in ncubes]#::Vector{typeof(mdbm.axes[i].ticks).parameters[1]}
for i in 1:length(mdbm.axes)]
end

"""
    getinterpolatedsolution(nc::NCube{IT,FT,N}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})

Provide the interpolated coordinates of the solution inside the provided n-cube `nc` (approximately where foo(x,y) == 0 and c(x,y)>0).
"""
function getinterpolatedsolution(nc::NCube{IT,FT,N}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    [
        (typeof(mdbm.axes[i].ticks).parameters[1])((mdbm.axes[i].ticks[nc.corner[i]] * (1.0 - (nc.posinterp[i] + 1.0) / 2.0) +
        mdbm.axes[i].ticks[nc.corner[i] + 1] * ((nc.posinterp[i] + 1.0) / 2.0)))#::Vector{typeof(mdbm.axes[i].ticks).parameters[1]}
for i in 1:length(mdbm.axes)]
end

"""
    getevaluatedpoints(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})

Provide all the coordinates where the function is evaluated.
"""
function getevaluatedpoints(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    [[x.callargs[i] for x in mdbm.fc.fvalarg] for i in 1:N]
end

"""
    getevaluatedpoints(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})

Provide the functionvalues for all the coordinates where the function is evaluated.
(see: `getevaluatedpoints`)
"""
function getevaluatedfunctionvalues(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    [x.funval for x in mdbm.fc.fvalarg]
end

"""
    getevaluatedconstraintvalues(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})

Provide the constraints values for all the coordinates where the function is evaluated.
(see: `getevaluatedpoints`)
"""
function getevaluatedconstraintvalues(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    [x.cval for x in mdbm.fc.fvalarg]
end

# function generateneighbours(ncubes::Vector{NCube},mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where N  where Nf where Nc
function generateneighbours(ncubes::Vector{NCube{IT,FT,N}}, mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    #TODO: Let the user select the method
    #-------all faceing/cornering neightbours----------------
    # #neighbourind=1:2^length(mdbm.axes) #cornering neightbours also - unnecessary
    # neighbourind=2 .^(0:(length(mdbm.axes)-1)) .+ 1 #neighbour only on the side
    # T101=[-mdbm.T01[neighbourind]...,mdbm.T01[neighbourind]...]
    #
    # nc_neighbour = Array{typeof(mdbm.ncubes[1])}(undef,0)
    # NumofNCubes=length(ncubes)
    # for iT in 1:length(T101)
    #     append!(nc_neighbour,deepcopy(ncubes))
    #     for nci in ((1+NumofNCubes*(iT-1)):(NumofNCubes+NumofNCubes*(iT-1)))
    #         nc_neighbour[nci].corner[:]=nc_neighbour[nci].corner+T101[iT].*nc_neighbour[nci].size
    #     end
    # end
    #-------all faceing/cornering neightbours----------------

    #-------direcational neightbours----------------
    nc_neighbour = Array{typeof(mdbm.ncubes[1])}(undef, 0)
    Nface = 2^N
    indpos = [[mdbm.T01[i][dir] for i in 1:Nface] for dir in 1:N]
    indneg = [[!mdbm.T01[i][dir] for i in 1:Nface] for dir in 1:N]

    for nc in ncubes
        fcvals = getcornerval(nc, mdbm)
        for dir in 1:N
            # indpos=[mdbm.T01[i][dir] for i in 1:Nface]
            # indneg=[!mdbm.T01[i][dir] for i in 1:Nface]
            if issingchange(fcvals[indpos[dir]], Nf, Nc)
                push!(nc_neighbour, deepcopy(nc))
                nc_neighbour[end].corner[dir] += nc_neighbour[end].size[dir]
            end
            if issingchange(fcvals[indneg[dir]], Nf, Nc)
                push!(nc_neighbour, deepcopy(nc))
                nc_neighbour[end].corner[dir] -= nc_neighbour[end].size[dir]
            end
        end
    end
    #-------direcational neightbours----------------
    filter!(nc->!(any(nc.corner .< 1) || any((nc.corner + nc.size) .> [length.(mdbm.axes)...])), nc_neighbour)#remove the overhanging ncubes
    sort!(nc_neighbour; alg = QuickSort)
    unique!(nc_neighbour)
    return nc_neighbour
end

"""
    checkneighbour!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}; interpolationorder::Int=0, maxiteration::Int=0)

Check the face-neighbouring n-cubes to discover the lost (missing) part of the solutions.

It is possible, that the sub-manifold (e.g.: curve) of the solution is not completely detected (it is interrupted).
This way the missing part can be explored similarly to a continuation method (with minimal function evaluation).

It works only if the dimension of the solution object is larger the zero (the maifold is not a set of poins).

    # Arguments
    - `interpolationorder::Int=0`: interpolation order method of the neighbours checked
    - `maxiteration::Int=0~: the max number of steps in the 'continuation-like' exploring. If zero, then infinity steps are allowed
"""
function checkneighbour!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}; interpolationorder::Int = 0, maxiteration::Int = 0) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}#only for unite size cubes

    if isempty(mdbm.ncubes)
        println("There is no bracketing n-cubes to check!")
    else
        ncubes2check = mdbm.ncubes
        numberofiteration = 0;
        while !isempty(ncubes2check) && (maxiteration == 0 ? true : numberofiteration < maxiteration)
            numberofiteration += 1
            ncubes2check = generateneighbours(ncubes2check, mdbm)

            deleteat!(ncubes2check, is_sorted_in_sorted(ncubes2check, mdbm.ncubes))#delete the ones which is already presented
            _interpolate!(ncubes2check, mdbm, Val{interpolationorder})#remove the non-bracketing, only proper new bracketing cubes remained

            append!(mdbm.ncubes, deepcopy(ncubes2check))
            sort!(mdbm.ncubes; alg = QuickSort)
        end
    end
end

"""
    connect(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT})

Provide the edge connection (as a list of point-paris) of the solution point-cloud based on the face-neighbouring n-cubes.
"""
function connect(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    #---------- line connection (no corener neighbour is needed) --------------
    DT1 = Array{Tuple{Int64,Int64}}(undef, 0)
    for inc in 1:length(mdbm.ncubes)
        ncneigh = generateneighbours([mdbm.ncubes[inc]], mdbm)
        indinresults = index_sorted_in_sorted(ncneigh, mdbm.ncubes)#delete the ones which is already presented
        append!(DT1, [(inc, x) for x  in indinresults if x != 0])
    end
    filter!(d->(d[2] > d[1]), DT1)#TODO: eleve csak ez egyik irányban levő szomszédokat kellene keresni!!!
    sort!(DT1; alg = QuickSort)
    unique!(DT1)
    return DT1
end

"""
    triangulation(DT1::Array{Tuple{Int64,Int64}})::Array{Tuple{Int64,Int64,Int64}}

Provide the triangulation (as a list of point-triplets) of the solution point-cloud based on the face-neighbouring n-cubes.
DT1 is the result of the edge connection performed by `connect`.
"""
function triangulation(DT1::Array{Tuple{Int64,Int64}})::Array{Tuple{Int64,Int64,Int64}}

    #DT=sort!()[DT1;[(d[2],d[1]) for d in DT1]]);
    DT = sort!([DT1;[d[[2,1]] for d in DT1]])#both direction of line connection is necessay!

    #L=[filter(d->d[1]==i,DT) for i in 1:length(mdbm.ncubes)]
    L = [
    [dd[2] for dd in filter(d->d[1] == i, DT)] for i in 1:maximum(DT1)[2]]

    DT4 = Array{Tuple{Int64,Int64,Int64,Int64}}(undef, 0)#quadratic patch
    for i in 1:size(L, 1)
        for j in L[i]#filter((L{i}>i) for i in )
            if i > j#i must be the larges value to remove the repetition of the some surface segment
                for k in L[j]#(L{j}>i) it seems to be slower
                    if i > k# there is no backstep, and i must be the largest   (Equivalent: %((k~=i) && (i>k))%back step is not allowed
                        for m in L[k]
                            if ((m != j) && (j > m)) #&& (i>m) #back step is not allowed, and i must be the largest
                                if any(i .== L[m])
                                    push!(DT4, (i, j, k, m))
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    return [[dt[[1,2,3]] for dt in DT4];[dt[[3,1,4]] for dt in DT4]]#DT2 triangular patch from the quadratic patch

end

"""
    solve!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, iteration::Int; interpolationorder::Int=1)

Refine the `MDBM_Problem` `iteration` times, then perform a neighbour check.
`interpolationorder` defines the interpolation method within the n-cubes.

# Examples
```julia
julia> solve!(mymdbm,4)
```
"""
function solve!(mdbm::MDBM_Problem{fcT,N,Nf,Nc,t01T,t11T,IT,FT,aT}, iteration::Int; interpolationorder::Int = 1) where {fcT ,N ,Nf ,Nc ,t01T ,t11T ,IT ,FT,aT}
    interpolate!(mdbm, interpolationorder = interpolationorder)
    for k = 1:iteration
        refine!(mdbm)
        interpolate!(mdbm, interpolationorder = interpolationorder)
    end
    checkneighbour!(mdbm)
    interpolate!(mdbm, interpolationorder = interpolationorder)
    return mdbm
end

function index_sorted_in_sorted(a::AbstractVector, b::AbstractVector)::Array{Int64,1}
    # index of a[i] in b (0 if not present)
    # a and b must be sorted
    containingindex = zeros(Int64, length(a))
    startindex = 1;
    for ind2check in 1:length(a)
        detectedrange = searchsorted(b[startindex:end], a[ind2check])
        startindex = max(detectedrange.stop, detectedrange.start) + startindex - 1
        containingindex[ind2check] = isempty(detectedrange) ? 0 : startindex
        if startindex > length(b)
            break
        end
    end
    return containingindex
end

function is_sorted_in_sorted(a::AbstractVector, b::AbstractVector)::Array{Bool,1}
    # is a[i] in b
    # a and b must be sorted
    iscontained = falses(length(a))
    startindex = 1;
    for ind2check in 1:length(a)
        detectedrange = searchsorted(b[startindex:end], a[ind2check])
        etectedrange = searchsorted(view(b,startindex:length(b)), a[ind2check]) # still: it is the critical line!!!
        iscontained[ind2check] = !isempty(detectedrange)
        startindex = max(detectedrange.stop, detectedrange.start) + startindex - 1
        if startindex > length(b)
            break
        end
    end
    return iscontained
end


text:
import NetworkViz
using LightGraphs
using ThreeJS
using Colors

export layout_spring, find_edges, drawWheel, drawGraph, drawGraphwithText, addEdge, removeEdge, addNode, removeNode, layout_spring, plot

function find_edges{T}(loc_x::Array{Float64,1},loc_y::Array{Float64,1},loc_z::Array{Float64,1},adj_matrix::Array{T,2})

    size(adj_matrix, 1) != size(adj_matrix, 2) && error("Adj. matrix must be square.")
    vertices = Tuple{Float64, Float64, Float64}[]
    const N = length(loc_x)
    for i = 1:N
        for j = 1:N
            i == j && continue
            if adj_matrix[i,j] != zero(eltype(adj_matrix))
                push!(vertices, (loc_x[i],loc_y[i],loc_z[i]), (loc_x[j],loc_y[j],loc_z[j]))
            end
        end
    end
    return vertices
end

function find_edges{T}(loc_x::Array{Float64,1},loc_y::Array{Float64,1},adj_matrix::Array{T,2})

    size(adj_matrix, 1) != size(adj_matrix, 2) && error("Adj. matrix must be square.")
    vertices = Tuple{Float64, Float64, Float64}[]
    const N = length(loc_x)
    for i = 1:N
        for j = 1:N
            i == j && continue
            if adj_matrix[i,j] != zero(eltype(adj_matrix))
                push!(vertices, (loc_x[i],loc_y[i],0), (loc_x[j],loc_y[j],0))
            end
        end
    end
    return vertices
end

function drawWheel(num::Int,t=1)
    g = WheelGraph(num)
    drawGraph(g,z=t)
end

function drawGraph(g::Union{LightGraphs.DiGraph,LightGraphs.Graph};
                               node=NodeProperty(Color[parse(Colorant,"#00004d") for i in 1:nv(g)],0.5,1),
                               edge=EdgeProperty("#ff3333",2),
                               z=1
                              )
    am = full(adjacency_matrix(g))
    loc_x, loc_y, loc_z = layout_spring(am,z)
    pts = zip(loc_x,loc_y,loc_z,node.color)
    if z == 1
        vertices = find_edges(loc_x, loc_y, loc_z, am)
    else
        vertices = find_edges(loc_x, loc_y, am)
    end
    plot(collect(pts),vertices, node, edge)
end

function drawGraph(am::Array{Int,2};
                            node=NodeProperty(Color[parse(Colorant,"#00004d") for i in 1:size(am,1)],0.5,1),
                            edge=EdgeProperty("#ff3333",2),
                            z=1
                           )
    loc_x, loc_y, loc_z = layout_spring(am,z)
    pts = zip(loc_x,loc_y,loc_z,node.color)
    if z == 1
        vertices = find_edges(loc_x, loc_y, loc_z, am)
    else
        vertices = find_edges(loc_x, loc_y, am)
    end
    plot(collect(pts),vertices, node, edge)
end

function addEdge(g::Graph, node1::Int, node2::Int, t=1)
    add_edge!(g,node1,node2)
    drawGraph(g,z=t)
end

function removeEdge(g::Graph, node1::Int, node2::Int, t=1)
    rem_edge!(g,node1,node2)
    drawGraph(g,z=t)
end

function addNode(g::Graph, t=1)
    add_vertex!(g)
    drawGraph(g,z=t)
end

function removeNode(g::Graph, node::Int, t=1)
    rem_vertex!(g,node)
    drawGraph(g,z=t)
end

function drawnode{T}(pts::Array{T,1},node::NodeProperty)
  if node.shape == 1
      ThreeJS.pointcloud(collect(pts)) <<
      [
        ThreeJS.pointmaterial(Dict(
        :color=>"white",
        :size=>node.size,
        :colorkind=>"vertex",
        :transparent=>true,
        :alphatest=>0.5,
        :texture=>"/pkg/NetworkViz/disc.png",
        ))
      ]
  else
      ThreeJS.pointcloud(collect(pts)) <<
      [
        ThreeJS.pointmaterial(Dict(
        :color=>"white",
        :size=>node.size,
        :colorkind=>"vertex",
        :transparent=>true,
        :alphatest=>0.5,
        ))
      ]
  end
end

function plot{T}(pts::Array{T,1}, vertices::Array{Tuple{Float64,Float64,Float64},1}, node::NodeProperty, edge::EdgeProperty)
  outerdiv() <<
  (
  initscene() <<
  [
      drawnode(pts,node),
      ThreeJS.line(vertices,kind="pieces") <<
      [
          ThreeJS.linematerial(Dict(
          :kind=>"basic",
          :color=>"$(edge.color)",
          :linewidth=>"$(edge.width)"
          ))
      ],
      pointlight(3.0, 3.0, 3.0),
      camera(0.0, 0.0, 5.0)
  ]
  )
end

function drawGraphwithText(g::Union{LightGraphs.DiGraph,LightGraphs.Graph},z=1)
    am = full(adjacency_matrix(g))
    loc_x, loc_y, loc_z = layout_spring(am,z)
    pts = zip(loc_x,loc_y,loc_z)
    if z == 1
        vertices = find_edges(loc_x, loc_y, loc_z, am)
    else
        vertices = find_edges(loc_x, loc_y, am)
    end
    outerdiv() <<
    (
    initscene() <<
    [
        ThreeJS.pointcloud(collect(pts)) <<
        [
          ThreeJS.pointmaterial(Dict(
          :color=>"#00004d",
          :size=>0.5,
          :transparent=>true,
          :alphatest=>0.2,
          :texture=>"/pkg/NetworkViz/disc.png",
          ))
        ],
        ThreeJS.line(vertices,kind="pieces") <<
        [
            ThreeJS.linematerial(Dict(
            :kind=>"basic",
            :color=>"#ff3333",
            :linewidth=>2
            ))
        ],
        [ThreeJS.text(map(x->x+0.08,item)...,"$idx") for (idx,item) in enumerate(pts)],
        pointlight(3.0, 3.0, 3.0),
        camera(0.0, 0.0, 5.0)
    ]
    )
end


text:
module Example

using Glimmer, Glimmer.FlexUI
using OpticSim, OpticSim.Geometry, OpticSim.Emitters
using OpticSimVis
using DataFrames
using StaticArrays
using TableView

println("Start [$(splitext(basename(@__FILE__))[1])]")

#---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
#---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Luxor")
prop!(app, :winInitWidth, 1400)
prop!(app, :winInitHeight,1200)

#---------------------------------------------------------------
# define the 3D Viewer
#---------------------------------------------------------------
scene = Scene(openWindow = false)
set_Z_up!(scene)
grid!(scene, false)
axes!(scene, false)
cameraTransform!(scene, lookAt(SVector(1.0, 1.0, 0.0)*50, SVector(0.0, 0.0, 0.0)*50))
cameraPlanes!(scene, 0.1, 1000.0)

#---------------------------------------------------------------
# Define Variables
#---------------------------------------------------------------
perscription_html = addVariable!(app, Variable(name="perscription_html", type="string",value=""))
radius1 = addVariable!(app, Variable(name="radius1", type="flota64",value=26))
radius2 = addVariable!(app, Variable(name="radius2", type="flota64",value=66))
count = addVariable!(app, Variable(name="count", type="flota64",value=16))

image = addVariable!(app, Variable(name="image", type="image", value="", ))

#---------------------------------------------------------------
# Define Controls
#---------------------------------------------------------------
ui = VContainer(
    Card(
        title="Controls",
        content=HContainer(
            Card(
                title="Radius 1",
                content=VContainer(
                    Slider(
                        text="Radius 1",
                        trailingText="[\$()]",
                        min=15,
                        max=40,
                        value=26,
                        variable="radius1"
                    ),
                    Field(
                        input="number",
                        label="Radius 1",
                        hint ="Save Field as the Slider above",
                        variable="radius1",
                    ),
                ),
            ),
            Card(
                title="Radius 2",
                content=VContainer(
                    Slider(
                        text="Radius 2",
                        trailingText="[\$()]",
                        min=15,
                        max=100,
                        value=26,
                        variable="radius2"
                    ),
                    Field(
                        input="number",
                        label="Radius 2",
                        hint ="Save Field as the Slider above",
                        variable="radius2",
                    ),
                ),
            ),
        ),
    ),

    Card(
        title="Perscription",
        # subtitle="Currently, Glimmer does not contain a GRID component, but allow you to utilize existing GRID component such as the one in the TableView package.",
        content=VContainer(
            RawHTML(
                html="\$(perscription_html)",
                style="width: 100%; height: 300px;"
            ),
        ),
    ),

    MeshCatViewer(
        url =  OpticSimVis.url(scene),
        width = "100%",
        height = "600px",
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control
)
# set the controls for the application
controls!(app, ui)

#---------------------------------------------------------------
# the render function - preparing the image
#---------------------------------------------------------------
function updatePerscription(df::DataFrame)
    table_data = TableView.showtable(df)
    perscription_html[] = renderHTML(table_data)
    Glimmer.FlexUI.forceUpdateControls!(app)
end

function render()
    clear(scene)

    Air = OpticSim.Air
    g1, g2 = OpticSim.SCHOTT.N_SK16, OpticSim.SCHOTT.N_SF2

    df = DataFrame(
        SurfaceType  = ["Object", "Standard", "Standard", "Standard", "Stop", "Standard", "Standard", "Image"],
        Radius       = [Inf,      26.777,     66.604,     -35.571,    35.571, 35.571,     -26.777,    Inf    ],
        Thickness    = [Inf,      4.0,        2.0,        4.0,        2.0,    4.0,        44.748,     missing],
        Material     = [Air,      g1,         Air,        g2,         Air,    g1,         Air,        missing],
        SemiDiameter = [Inf,      8.580,      7.513,      7.054,      6.033,  7.003,      7.506,      15.0   ],
    )
    df[2,"Radius"] = radius1[]
    df[3,"Radius"] = radius2[]

    updatePerscription(df)

    sys = AxisymmetricOpticalSystem{Float64}(df)

    origins = Origins.Hexapolar(8, 15.0, 15.0)
    directions = Directions.Constant(0.0, 0.0, -1.0)
    s1 = Sources.Source(; origins, directions, sourcenum=1)

    transform = Transform(rotmatd(10, 0, 0), Geometry.unitZ3())
    s2 = Sources.Source(; transform, origins, directions, sourcenum=2)

    raygenerator = Sources.CompositeSource(Transform(), [s1, s2])

    trackallrays = test = colorbysourcenum = true; resolution = (1000, 700)

    # @info "-"^50
    # @info typeof(sys)
    # @info typeof(sys.system)
    # @info typeof(sys.system.assembly)
    # @info "-"^50
    OpticSimVis.draw!(scene, sys.system)

    resetdetector!(sys)
    if (true)
        OpticSimVis.drawtracerays!(
            scene, sys;
            raygenerator=raygenerator,
            trackallrays = true,
            colorbynhits = true,
            test = true,
            numdivisions = 100,
            drawgen = false,
            drawsys = true
        )
    end

    # Vis.drawtracerays(sys; raygenerator, trackallrays, test, colorbysourcenum, resolution)
    # Vis.make2dy(); Vis.save(filename)
    # return sys

end
renderFunction!(app, render)

#---------------------------------------------------------------
# Run the application
#---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module


text:
struct HH3DDoubleLayerNear{K}
  gamma::K
end

function HH3DDoubleLayerNear(;wavenumber=error("wavenumber is a required argument"))
  if iszero(real(wavenumber))
    HH3DDoubleLayerNear(-imag(wavenumber))
  else
    HH3DDoubleLayerNear(wavenumber*im)
  end
end

quaddata(op::HH3DDoubleLayerNear,rs,els) = quadpoints(rs,els,(3,))
quadrule(op::HH3DDoubleLayerNear,refspace,p,y,q,el,qdata) = qdata[1,q]

function kernelvals(op::HH3DDoubleLayerNear,y,p)

  γ = op.gamma
  x = cartesian(p)
  r = y - x
  R = norm(r)
  γR = γ*R

  inv_R = 1/R

  expn = exp(-γR)
  green = expn * inv_R * inv_4pi
  gradgreen = -(γ + inv_R) * green * inv_R * r

  nx = normal(p)

  (;γ, r, R, green, gradgreen, nx)
end

function integrand(op::HH3DDoubleLayerNear,krn,y,f,p)

  ∇G = krn.gradgreen
  nx = krn.nx

  fx = f.value

  ∂G∂n = nx ⋅ ∇G

  return ∂G∂n * fx
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 57aa3170-6e51-11eb-2670-9b77065bc4da
begin
    using PlutoUI
    using Plots
    using Printf
	using SpecialFunctions
end

# ╔═╡ ad3d9470-6e82-11eb-3072-43decfc60ef9
md"_MAC Method for Lid-driven Cavity Problem, version 1_"

# ╔═╡ 484fcd20-6f23-11eb-0dc8-5d49530983a1
md"""
## **Lid-driven Cavity Problem with MAC Method**
#### **Description:**
This notebook uses MAC method to solve incompressible 2D Navier-Stokes equations for the solution of the lid-driven cavity problem.

		MAC(Re,h)

If the arguments are missing, it will use the default values.

		MAC(Re=0.002,h=0.002)

#### **Lid-driven Cavity Problem**
In this problem, there is a moving lid over a cavity. The lid is moving with a constant speed. The governing equations for the system can be written in the conservative form:

$$\frac{∂u}{∂x}+\frac{∂u}{∂y}=0$$
$$\frac{∂u}{∂t}+\frac{∂}{∂x}(u^2+\frac{p}{ρ})+\frac{∂(uv)}{∂y}=ν(\frac{∂^2u}{∂x^2}+\frac{∂^2u}{∂y^2})$$
$$\frac{∂v}{∂t}+\frac{∂(uv)}{∂x}(v^2+\frac{p}{ρ})+\frac{∂(uv)}{∂y}=ν(\frac{∂^2v}{∂x^2}+\frac{∂^2v}{∂y^2})$$

These equations can be nondimensionalized with the following values:

$$Re=\frac{ρ_\infty u_\infty L}{μ_\infty}$$
$$t^*=\frac{t u_\infty}{L}$$
$$x^*=\frac{x}{L}$$
$$y^*=\frac{y}{L}$$
$$u^*=\frac{u}{u_\infty}$$
$$v^*=\frac{v}{u_\infty}$$
$$p^*=\frac{p}{ρ_\infty u^2_\infty}$$

After the nondimensionalization * is ignored for simplicity:

$$\frac{∂u}{∂x}+\frac{∂u}{∂y}=0$$
$$\frac{∂u}{∂t}+\frac{∂}{∂x}(u^2+p)+\frac{∂(uv)}{∂y}=\frac{1}{Re}(\frac{∂^2u}{∂x^2}+\frac{∂^2u}{∂y^2})$$
$$\frac{∂v}{∂t}+\frac{∂(uv)}{∂x}(v^2+p)+\frac{∂(uv)}{∂y}=\frac{1}{Re}(\frac{∂^2v}{∂x^2}+\frac{∂^2v}{∂y^2})$$

Since continuity term does not include a time dependent term, there is not communication mechanism between continuity and momentum equations. To overcome this problem, possion equation will be solved. In this notebook, calculations are done in three steps:

1) Calculation of u*.
2) Calculation of p.
3) Calculation of corrected u.

This approach is called pressure correction. The discritization for the system is:

$$I   = \frac{(u_{i+1,j}+u_{i,j})^2-(u_{i,j}+u_{i-1,j})^2}{4Δx}$$
$$II  = \frac{(u_{i,j+1}+u_{i,j})(v_{i+1,j}+v_{i,j})}{4Δy}-\frac{(u_{i,j}+u_{i,j-1})(v_{i+1,j-1}+v_{i,j-1})}{4Δy}$$
$$III = \frac{1}{Re}(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{Δx^2}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{Δy^2})$$
$$\frac{u_{i,j}^*-u_{i,j}^n}{Δt}=-I-II+III$$

For the pressure correction:

$$I=\frac{p_{i+1,j}-2p_{i,j}+p_{i-1,j}}{Δx^2}+\frac{p_{i,j+1}-2p_{i,j}+p_{i,j-1}}{Δy^2}$$
$$II=\frac{\frac{u_{i,j}^*-u_{i-1,j}^*}{Δx}+\frac{v_{i,j}^*-v_{i,j-1}^*}{Δy}}{Δt}$$
$$p_{i,j}=p_{i,j}+ω\frac{Δx^2+Δy^2}{2}(I-II)$$

In this notebook, ω is taken as 1 and this correction method is applied 10 times. Once the pressure is calculated, correct velocity can be calculated for the next time step as:

$$\frac{u^{n+1}_{i,j}-u^*_{i,j}}{Δt}=\frac{p_{i+1,j}-p_{i,j}}{Δx}$$

Same procedure can be applied for the v velocity. Boundary conditions for the system is:

$$u(0,:)=0$$
$$u(N,:)=0$$
$$u(:,0)=0$$
$$u(:,N)=1$$
$$v(0,:)=0$$
$$v(N,:)=0$$
$$v(:,0)=0$$
$$v(:,N)=0$$

For the pressure Nuemann boundary condition is applied as:

$$\frac{∂p}{∂η}=0$$

where η is the wall normal direction.

Details of schemes:
Computational Fluid Dynamics, Volume I, Hoffmann - Chiang

Feel free to ask questions!

*Furkan Oz*
foz@okstate.edu
"""

# ╔═╡ 6bc808d0-6e51-11eb-0794-b9585413d41b

"""
This notebook uses MAC method to solve incompressible 2D Navier-Stokes equations for the solution of the lid-driven cavity problem.

		MAC(Re,h)

If the arguments are missing, it will use the default values.

		MAC(Re=0.002,h=0.002)

Furkan Oz,
foz@okstate.edu,

"""
function MAC(Re=100, h=0.02)
	# Initialiazing for domain setup
	xdom = [0 1.0] # Domain limits for x
	ydom = [0 1.0] # Domain limits for y
	Δx = h
	Δy = h
	Δx²= Δx^2
	Δy²= Δy^2
	nex = floor(Int32,(xdom[2]-xdom[1])/h) #Number of elements in x-direction
	ney = floor(Int32,(xdom[2]-xdom[1])/h) #Number of elements in y-direction
	nnx = nex + 1 #Number of nodes in x-direction
	nny = ney + 1 #Number of nodes in y-direction

	# Initialiazing for solver setup

	CFL = 1 # CFL condition (less than 1 preferred)
	Δtc = CFL*min(Δx,Δy) # Time-step based on Convective Term
	Δtv = 0.25*Re*min(Δx²,Δy²) # Time-step based on Viscous Term
	σ = 0.5
	a =-2.0/(Δx²) - 2.0/(Δy²)
	Δt =  σ*min(Δtc,Δtv) # Time-step based on lowest condition
	ϵ = 1e-5 # Convergence limit

	ttmax = 30000 # Maximum iteration limit
	iimax = 10 # Pressure correction iteration limit

	u  = zeros(nnx,nny+1)
	v  = zeros(nnx+1,nny)
	uˢ = zeros(nnx,nny+1)
	vˢ = zeros(nnx+1,nny)
	u₀ = zeros(nnx,nny+1)
	v₀ = zeros(nnx+1,nny)

	p  = zeros(nex+2,ney+2)
	p₀ = zeros(nex+2,ney+2)

	rᵤ = 0
	rᵥ = 0
	rₚ = 0
	ttend = 0

	rᵤₜ = zeros(ttmax)
	rᵥₜ = zeros(ttmax)
	rₚₜ = zeros(ttmax)

	#Solver starts

	for tt=1:ttmax
		for j=1:nny+1
			for i=1:nnx
				u₀[i,j] = u[i,j]
			end
		end
		for j=1:nny
			for i=1:nnx+1
				v₀[i,j] = v[i,j]
			end
		end
		for j=1:ney+2
			for i=1:nex+2
				p₀[i,j] = p[i,j]
			end
		end

		# Boundary conditions for velocity
		for i=1:nnx
			u[i,1] =-u[i,2]
			u[i,nny+1] = 2.0-u[i,nny]
		end
		for i=1:nny
			v[1,i] =-v[2,i]
			v[nnx+1,i] =-v[nnx,i]
		end

		# uˢ calculation
		for j=2:nny
			for i=2:nnx-1

  				c₃ = (1.0/Re)*((u[i+1,j]-2.0*u[i,j]+u[i-1,j])/(Δx²)  				                      + (u[i,j+1]-2.0*u[i,j]+u[i,j-1])/(Δy²) )

  				c₁ = (0.25/Δx)*((u[i+1,j]+u[i,j])^2 - (u[i,j]+u[i-1,j])^2)

  				c₂ = (0.25/Δy)*((u[i,j+1]+u[i,j])*(v[i+1,j]+v[i,j])                                      - (u[i,j]+u[i,j-1])*(v[i+1,j-1]+v[i,j-1]))

  				uˢ[i,j] = u[i,j] + Δt*(-c₁-c₂+c₃)

			end
		end

		# vˢ calculation
		for j=2:nny-1
			for i=2:nnx

  				c₃ = (1.0/Re)*((v[i+1,j]-2.0*v[i,j]+v[i-1,j])/(Δx²)  				                      + (v[i,j+1]-2.0*v[i,j]+v[i,j-1])/(Δy²) )

  				c₂ = (0.25/Δy)*((v[i,j+1]+v[i,j])^2 - (v[i,j]+v[i,j-1])^2)

  				c₁ = (0.25/Δx)*((u[i,j+1]+u[i,j])*(v[i+1,j]+v[i,j])                                      - (u[i-1,j+1]+u[i-1,j])*(v[i,j]+v[i-1,j]))

  				vˢ[i,j] = v[i,j] + Δt*(-c₁-c₂+c₃)

			end
		end

		# Boundary conditions for velocity correction
		for i=1:nnx
			uˢ[i,1] =-uˢ[i,2]
			uˢ[i,nny+1] = 2.0-uˢ[i,nny]
		end
		for i=1:nny
			vˢ[1,i] =-vˢ[2,i]
			vˢ[nnx+1,i] =-vˢ[nnx,i]
		end

		# Boundary conditions for pressure
		for i=1:nex+2
			p[i,1] = p[i,2]
			p[i,ney+2] = p[i,ney+1]
		end
		for i=1:ney+2
			p[1,i] = p[2,i]
			p[nex+2,i] = p[nex+1,i]
		end

		# Poisson Solver
		ω = 1.0
		for ii = 1:iimax
			for j=2:ney+1
				for i=2:nex+1

  					f = ((uˢ[i,j]-uˢ[i-1,j])/Δx + (vˢ[i,j]-vˢ[i,j-1])/Δy)/Δt

  					r = (f - (p[i+1,j]-2.0*p[i,j]+p[i-1,j])/(Δx²)                                           -      (p[i,j+1]-2.0*p[i,j]+p[i,j-1])/(Δy²) )

  					p[i,j] = p[i,j] + ω*r/a

				end
			end
		end

		# Correct velocity calculation
		for j=2:nny
			for i=2:nnx-1
  				u[i,j] = uˢ[i,j] - Δt*(p[i+1,j]-p[i,j])/Δx
			end
		end

		for j=2:nny-1
			for i=2:nnx
  				v[i,j] = vˢ[i,j] - Δt*(p[i,j+1]-p[i,j])/Δy
			end
		end

		#check residual
		rᵤ = 0.0
		rᵥ = 0.0
		rₚ = 0.0
		for j=2:nny
			for i=2:nnx
  				rᵤ = rᵤ + (u[i,j] - u₀[i,j])^2
  				rᵥ = rᵥ + (v[i,j] - v₀[i,j])^2
  				rₚ = rₚ + (p[i,j] - p₀[i,j])^2
			end
		end
		rᵤ = sqrt(rᵤ/(nex*ney))
		rᵥ = sqrt(rᵥ/(nex*ney))
		rₚ = sqrt(rₚ/(nex*ney))

		rᵤₜ[tt] = rᵤ
		rᵥₜ[tt] = rᵥ
		rₚₜ[tt] = rₚ
		ttend = tt
		@show tt rᵤ rᵥ rₚ ttend rᵤₜ[tt]
		if (rᵤ<=ϵ && rᵥ<=ϵ && rₚ<=ϵ)
			break
		end

	end
	ttend = 1:ttend

	return u,v,p,nnx,ttend,rᵤₜ,rᵥₜ,rₚₜ
end

# ╔═╡ 9be09b70-6f2a-11eb-0f64-018b8c98e73f
begin
	u100c,v100c,p100c = MAC(100)
	contourf(u100c',levels=50,c= :redsblues, clims=(-1, 1),linewidth = 0,title = "Re=100 - u Contour")
end

# ╔═╡ 6c9f9cc0-6e54-11eb-380f-9f8f2d11343c
begin
	Re100data = [1	1	0.9766	0.84123	0.9688	0.78871	0.9609	0.73722	0.9531	0.68717	0.8516	0.23151	0.7344	0.00332	0.6172	-0.13641	0.5000 -0.20581	0.4531 -0.21090	0.2813	-0.15662	0.1719	-0.1015	0.1016	-0.06434	0.0703	-0.04775	0.0625	-0.04192	0.0547	-0.03717	0	0]
	y100data = [Re100data[i] for i=1:2:34]
	u100data = [Re100data[i] for i=2:2:34]

	u100_0_02,v100_0_02,p100_0_02,nnx100_0_02 = MAC(100,0.02)
	y100_0_02 = -0.01:0.02:1+0.01
	nnx100_0_02m = floor(Int32,(nnx100_0_02-1)/2)

	u100_0_01,v100_0_01,p100_0_01,nnx100_0_01 = MAC(100,0.01)
	y100_0_01 = -0.005:0.01:1+0.005
	nnx100_0_01m = floor(Int32,(nnx100_0_01-1)/2)

	u100_0_005,v100_0_005,p100_0_005,nnx100_0_005 = MAC(100,0.005)
	y100_0_005 = -0.0025:0.005:1+0.0025
	nnx100_0_005m = floor(Int32,(nnx100_0_005-1)/2)

	plot([u100_0_02[nnx100_0_02m,:], u100_0_01[nnx100_0_01m,:], u100_0_005[nnx100_0_005m,:], u100data],[y100_0_02,y100_0_01,y100_0_005,y100data],
				title = "Re=100 - Center u Velocity",
				label = ["MAC h=0.02" "MAC h=0.01" "MAC h=0.005" "Ghia"],
				legend = :bottomright,
				#ylims = (0,0.04),
				xlabel = "U",
				ylabel = "y",
				linewidth = 2,
				linecolor = [:black :black :black :transparent],
				markershape = :auto,
				markercolor = [:red :red :red :green],
				)
end

# ╔═╡ 07e46510-6f1f-11eb-21a0-f50e1e1e67b8
begin

	 x,x,x,x,tt100_0_02,ru100_0_02,rv100_0_02,rp100_0_02 = MAC(100,0.02)

	plot([tt100_0_02],[log.(ru100_0_02[1:tt100_0_02[end]]),log.(rv100_0_02[1:tt100_0_02[end]]),log.(rp100_0_02[1:tt100_0_02[end]])],
				title = "Re=100 - h=0.02 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				xlims = (0,1500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)
end

# ╔═╡ a71f5160-6f26-11eb-1ffb-15467ef1e901
begin

	q,q,q,q,tt100_0_01,ru100_0_01,rv100_0_01,rp100_0_01 = MAC(100,0.01)

	plot([tt100_0_01],[log.(ru100_0_01[1:tt100_0_01[end]]),log.(rv100_0_01[1:tt100_0_01[end]]),log.(rp100_0_01[1:tt100_0_01[end]])],
				title = "Re=100 - h=0.01 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)
end

# ╔═╡ 13382bd0-6f2a-11eb-392a-6fc93389bda6
begin

	w,w,w,w,tt100_0_005,ru100_0_005,rv100_0_005,rp100_0_005 = MAC(100,0.005)

	plot([tt100_0_005],[log.(ru100_0_005[1:tt100_0_005[end]]),log.(rv100_0_005[1:tt100_0_005[end]]),log.(rp100_0_005[1:tt100_0_005[end]])],
				title = "Re=100 - h=0.005 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)
end

# ╔═╡ 4c463b00-6f2b-11eb-2840-7384b467a27b
begin
	u400c,v400c,p400c = MAC(400)
	contourf(u400c',levels=50,c= :redsblues, clims=(-1, 1),linewidth = 0,title = "Re=400 - u Contour")
end

# ╔═╡ bfad9c20-6e7f-11eb-2fc9-49f1df970d96
begin
	Re400data = [1	1	0.9766	0.75837	0.9688	0.68439	0.9609	0.61756	0.9531	0.55892	0.8516	0.29093	0.7344	0.16256	0.6172	0.02135	0.5000 -0.11477	0.4531 -0.17119	0.2813	-0.32726	0.1719	-0.24299	0.1016	-0.14612	0.0703	-0.10338	0.0625	-0.09266	0.0547	-0.08186	0	0
]
	y400data = [Re400data[i] for i=1:2:34]
	u400data = [Re400data[i] for i=2:2:34]

	u400_0_02,v400_0_02,p400_0_02,nnx400_0_02 = MAC(400,0.02)
	y400_0_02 = -0.01:0.02:1+0.01
	nnx400_0_02m = floor(Int32,(nnx400_0_02-1)/2)

	u400_0_01,v400_0_01,p400_0_01,nnx400_0_01 = MAC(400,0.01)
	y400_0_01 = -0.005:0.01:1+0.005
	nnx400_0_01m = floor(Int32,(nnx400_0_01-1)/2)

	u400_0_005,v400_0_005,p400_0_005,nnx400_0_005 = MAC(400,0.005)
	y400_0_005 = -0.0025:0.005:1+0.0025
	nnx400_0_005m = floor(Int32,(nnx400_0_005-1)/2)

	plot([u400_0_02[nnx400_0_02m,:], u400_0_01[nnx400_0_01m,:], u400_0_005[nnx400_0_005m,:], u400data],[y400_0_02,y400_0_01,y400_0_005,y400data],
				title = "Re=400 - Center Velocity",
				label = ["MAC h=0.02" "MAC h=0.01" "MAC h=0.005" "Ghia"],
				legend = :bottomright,
				#ylims = (0,0.04),
				xlabel = "U",
				ylabel = "y",
				linewidth = 2,
				linecolor = [:black :black :black :transparent],
				markershape = :auto,
				markercolor = [:red :red :red :green],
				)
end

# ╔═╡ 26043960-6f2b-11eb-2fdd-33f7de7fc9e3
begin

	 t,t,t,t,tt400_0_02,ru400_0_02,rv400_0_02,rp400_0_02 = MAC(400,0.02)

	plot([tt400_0_02],[log.(ru400_0_02[1:tt400_0_02[end]]),log.(rv400_0_02[1:tt400_0_02[end]]),log.(rp400_0_02[1:tt400_0_02[end]])],
				title = "Re=400 - h=0.02 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				xlims = (0,1500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)
end

# ╔═╡ 03cf5b90-6f2b-11eb-0dd8-7f2061f9d9d7
begin

	r,r,r,r,tt400_0_01,ru400_0_01,rv400_0_01,rp400_0_01 = MAC(400,0.01)

	plot([tt400_0_01],[log.(ru400_0_01[1:tt400_0_01[end]]),log.(rv400_0_01[1:tt400_0_01[end]]),log.(rp400_0_01[1:tt400_0_01[end]])],
				title = "Re=400 - h=0.01 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)
end

# ╔═╡ e680ea90-6f2a-11eb-1ca2-89a396add38a
begin

	e,e,e,e,tt400_0_005,ru400_0_005,rv400_0_005,rp400_0_005 = MAC(400,0.005)

	plot([tt400_0_005],[log.(ru400_0_005[1:tt400_0_005[end]]),log.(rv400_0_005[1:tt400_0_005[end]]),log.(rp400_0_005[1:tt400_0_005[end]])],
				title = "Re=400 - h=0.005 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)
end

# ╔═╡ 57ed04c2-6f2b-11eb-20c8-1725c926dd2a
begin
	u1000c,v1000c,p1000c = MAC(1000)
	contourf(u1000c',levels=50,c= :redsblues, clims=(-1, 1),linewidth = 0,title = "Re=1000 - u Contour")
end

# ╔═╡ c8ece650-6e80-11eb-2882-7fed169d61ed
begin
	Re1000data = [1	1	0.9766	0.65928	0.9688	0.57492	0.9609	0.51117	0.9531	0.46604	0.8516	0.33304	0.7344	0.18719	0.6172	0.05702	0.5	-0.0608	0.4531	-0.10648	0.2813	-0.27805	0.1719	-0.38289	0.1016	-0.2973	0.0703	-0.2222	0.0625	-0.20196	0.0547	-0.18109	0	0]
	y1000data = [Re1000data[i] for i=1:2:34]
	u1000data = [Re1000data[i] for i=2:2:34]

	u1000_0_02,v1000_0_02,p1000_0_02,nnx1000_0_02 = MAC(1000,0.02)
	y1000_0_02 = -0.01:0.02:1+0.01
	nnx1000_0_02m = floor(Int32,(nnx1000_0_02-1)/2)

	u1000_0_01,v1000_0_01,p1000_0_01,nnx1000_0_01 = MAC(1000,0.01)
	y1000_0_01 = -0.005:0.01:1+0.005
	nnx1000_0_01m = floor(Int32,(nnx1000_0_01-1)/2)

	u1000_0_005,v1000_0_005,p1000_0_005,nnx1000_0_005 = MAC(1000,0.005)
	y1000_0_005 = -0.0025:0.005:1+0.0025
	nnx1000_0_005m = floor(Int32,(nnx1000_0_005-1)/2)

	plot([u1000_0_02[nnx1000_0_02m,:], u1000_0_01[nnx1000_0_01m,:], u1000_0_005[nnx1000_0_005m,:], u1000data],[y1000_0_02,y1000_0_01,y1000_0_005,y1000data],
				title = "Re=1000 - Center Velocity",
				label = ["MAC h=0.02" "MAC h=0.01" "MAC h=0.005" "Ghia"],
				legend = :bottomright,
				#ylims = (0,0.04),
				xlabel = "U",
				ylabel = "y",
				linewidth = 2,
				linecolor = [:black :black :black :transparent],
				markershape = :auto,
				markercolor = [:red :red :red :green],
				)
end

# ╔═╡ 94b982e0-6f29-11eb-2f68-839a418ad1f8
begin

	 c,c,c,c,tt1000_0_02,ru1000_0_02,rv1000_0_02,rp1000_0_02 = MAC(1000,0.02)

	plot([tt1000_0_02],[log.(ru1000_0_02[1:tt1000_0_02[end]]),log.(rv1000_0_02[1:tt1000_0_02[end]]),log.(rp1000_0_02[1:tt1000_0_02[end]])],
				title = "Re=1000 - h=0.02 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,1500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)

end

# ╔═╡ b5637060-6f28-11eb-0850-5dff9a108c4b
begin

	 b,b,b,b,tt1000_0_01,ru1000_0_01,rv1000_0_01,rp1000_0_01 = MAC(1000,0.01)

	plot([tt1000_0_01],[log.(ru1000_0_01[1:tt1000_0_01[end]]),log.(rv1000_0_01[1:tt1000_0_01[end]]),log.(rp1000_0_01[1:tt1000_0_01[end]])],
				title = "Re=1000 - h=0.01 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,1500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)

end

# ╔═╡ 2e0fecf0-6f29-11eb-1d52-71a1c891ed5a
begin

	 a,a,a,a,tt1000_0_005,ru1000_0_005,rv1000_0_005,rp1000_0_005 = MAC(1000,0.005)

	plot([tt1000_0_005],[log.(ru1000_0_005[1:tt1000_0_005[end]]),log.(rv1000_0_005[1:tt1000_0_005[end]]),log.(rp1000_0_005[1:tt1000_0_005[end]])],
				title = "Re=1000 - h=0.005 - Residual",
				label = ["u residual" "v residual" "p residual"],
				legend = :topright,
				#xlims = (0,1500),
				ylims = (-15,-2.5),
				xlabel = "Iteration",
				ylabel = "log(Residual)",
				linewidth = 2,
				linecolor = [:auto],
				#markershape = :auto,
				#markercolor = [:red],
				)

end

# ╔═╡ Cell order:
# ╟─57aa3170-6e51-11eb-2670-9b77065bc4da
# ╟─ad3d9470-6e82-11eb-3072-43decfc60ef9
# ╟─484fcd20-6f23-11eb-0dc8-5d49530983a1
# ╟─6bc808d0-6e51-11eb-0794-b9585413d41b
# ╟─9be09b70-6f2a-11eb-0f64-018b8c98e73f
# ╟─6c9f9cc0-6e54-11eb-380f-9f8f2d11343c
# ╟─07e46510-6f1f-11eb-21a0-f50e1e1e67b8
# ╟─a71f5160-6f26-11eb-1ffb-15467ef1e901
# ╟─13382bd0-6f2a-11eb-392a-6fc93389bda6
# ╟─4c463b00-6f2b-11eb-2840-7384b467a27b
# ╟─bfad9c20-6e7f-11eb-2fc9-49f1df970d96
# ╟─26043960-6f2b-11eb-2fdd-33f7de7fc9e3
# ╟─03cf5b90-6f2b-11eb-0dd8-7f2061f9d9d7
# ╟─e680ea90-6f2a-11eb-1ca2-89a396add38a
# ╟─57ed04c2-6f2b-11eb-20c8-1725c926dd2a
# ╟─c8ece650-6e80-11eb-2882-7fed169d61ed
# ╟─94b982e0-6f29-11eb-2f68-839a418ad1f8
# ╟─b5637060-6f28-11eb-0850-5dff9a108c4b
# ╟─2e0fecf0-6f29-11eb-1d52-71a1c891ed5a


text:
# Common functions

# KNITRO special types
const KNLONG = Clonglong
const KNBOOL = Cint

"A macro to make calling KNITRO's KN_* C API a little cleaner"
macro kn_ccall(func, args...)
    f = Base.Meta.quot(Symbol("KN_$(func)"))
    args = [esc(a) for a in args]
    quote
        ccall(($f, libknitro), $(args...))
    end
end

"Check if return value is valid."
function _checkraise(ret::Cint)
    if ret != 0
        error("Fail to use specified function: $ret")
    end
end

"Format output returned by KNITRO as proper Julia string."
function format_output(output::AbstractString)
    # remove trailing whitespace
    res = strip(output)
    # remove special characters
    res = strip(res, '\0')
end

"Return the current KNITRO version."
function get_release()
    len = 15
    out = zeros(Cchar,len)

    @kn_ccall(get_release, Cvoid, (Cint, Ptr{Cchar}), len, out)
    return String(strip(String(convert(Vector{UInt8},out)), '\0'))
end


text:
@enum DUCKDB_TYPE_::UInt32 begin
    DUCKDB_TYPE_INVALID = 0
    DUCKDB_TYPE_BOOLEAN
    DUCKDB_TYPE_TINYINT
    DUCKDB_TYPE_SMALLINT
    DUCKDB_TYPE_INTEGER
    DUCKDB_TYPE_BIGINT
    DUCKDB_TYPE_UTINYINT
    DUCKDB_TYPE_USMALLINT
    DUCKDB_TYPE_UINTEGER
    DUCKDB_TYPE_UBIGINT
    DUCKDB_TYPE_FLOAT
    DUCKDB_TYPE_DOUBLE
    DUCKDB_TYPE_TIMESTAMP
    DUCKDB_TYPE_DATE
    DUCKDB_TYPE_TIME
    DUCKDB_TYPE_INTERVAL
    DUCKDB_TYPE_HUGEINT
    DUCKDB_TYPE_VARCHAR
    DUCKDB_TYPE_BLOB
end

const DUCKDB_TYPE = DUCKDB_TYPE_

"""
Days are stored as days since 1970-01-01\n
Use the duckdb_from_date/duckdb_to_date function to extract individual information

"""
struct duckdb_date
    days::Int32
end

struct duckdb_date_struct
    year::Int32
    month::Int8
    day::Int8
end

"""
Time is stored as microseconds since 00:00:00\n
Use the duckdb_from_time/duckdb_to_time function to extract individual information

"""
struct duckdb_time
    micros::Int64
end

struct duckdb_time_struct
    hour::Int8
    min::Int8
    sec::Int8
    micros::Int32
end

"""
Timestamps are stored as microseconds since 1970-01-01\n
Use the duckdb_from_timestamp/duckdb_to_timestamp function to extract individual information

"""
struct duckdb_timestamp
    micros::Int64
end

struct duckdb_timestamp_struct
    date::Ref{duckdb_date_struct}
    time::Ref{duckdb_time_struct}
end

struct duckdb_interval
    months::Int32
    days::Int32
    micros::Int64
end

"""
Hugeints are composed in a (lower, upper) component\n
The value of the hugeint is upper * 2^64 + lower\n
For easy usage, the functions duckdb_hugeint_to_double/duckdb_double_to_hugeint are recommended

"""
struct duckdb_hugeint
    lower::UInt64
    upper::Int64
end

struct duckdb_blob
    data::Cvoid
    size::UInt64
end

struct duckdb_column
    data::Ptr{Cvoid}
    nullmask::Ptr{UInt8}
    type::Ptr{DUCKDB_TYPE}
    name::Ptr{UInt8}
    internal_data::Ptr{Cvoid}
end

struct duckdb_result
    column_count::Ptr{UInt64}
    row_count::Ptr{UInt64}
    rows_changed::Ptr{UInt64}
    columns::Ptr{duckdb_column}
    error_message::Ptr{UInt8}
    internal_data::Ptr{Cvoid}
end

DUCKDB_TYPES = Dict(
    DUCKDB_TYPE_BOOLEAN => UInt8,
    DUCKDB_TYPE_TINYINT => Int8,
    DUCKDB_TYPE_SMALLINT => Int16,
    DUCKDB_TYPE_INTEGER => Int32,
    DUCKDB_TYPE_BIGINT => Int64,
    DUCKDB_TYPE_UTINYINT => UInt8,
    DUCKDB_TYPE_USMALLINT => UInt16,
    DUCKDB_TYPE_UINTEGER => UInt32,
    DUCKDB_TYPE_UBIGINT => UInt64,
    DUCKDB_TYPE_FLOAT => Float32,
    DUCKDB_TYPE_DOUBLE => Float64,
    DUCKDB_TYPE_TIMESTAMP => Int64,
    DUCKDB_TYPE_DATE => Int32,
    DUCKDB_TYPE_TIME => Int64,
    DUCKDB_TYPE_INTERVAL => duckdb_interval,
    DUCKDB_TYPE_HUGEINT => duckdb_hugeint,
    DUCKDB_TYPE_VARCHAR => Ptr{UInt8},
    DUCKDB_TYPE_BLOB => duckdb_blob,
)


text:
module test_lupin_vecjoin

using Test
using Lupin

@test vecjoin([], +) == []
@test vecjoin([1], +) == [1]
@test vecjoin([1, 2, 3], +) == [1, +, 2, +, 3]
@test vecjoin([1 2 3], +) == [1, +, 2, +, 3]
@test vecjoin([1 2 3], 5) == [1, 5, 2, 5, 3]
@test vecjoin(1:3, +) == [1, +, 2, +, 3]

@test isequal(vecjoin([1, 2, 3], NaN), [1, NaN, 2, NaN, 3])

end # module test_lupin_vecjoin