Domains.jl

module Domains

import ..Geometry: Geometry, float_type
using IntervalSets
export RectangleDomain

"""
    AbstractDomain

A domain represents a region of space.
"""
abstract type AbstractDomain end

function Base.summary(io::IO, domain::AbstractDomain)
    print(io, nameof(typeof(domain)))
end

const BCTagType = Union{Nothing, Tuple{Symbol, Symbol}}

float_type(domain::AbstractDomain) = float_type(coordinate_type(domain))

"""
    boundary_names(obj::Union{AbstractDomain, AbstractMesh, AbstractTopology})

A tuple or vector of unique boundary names of a spatial domain.
"""
function boundary_names end

struct IntervalDomain{CT, B} <: AbstractDomain where {
    CT <: Geometry.Abstract1DPoint{FT},
    B <: BCTagType,
} where {FT}
    coord_min::CT
    coord_max::CT
    boundary_names::B
end

const XIntervalDomain{FT, B} = IntervalDomain{Geometry.XPoint{FT}, B}
const YIntervalDomain{FT, B} = IntervalDomain{Geometry.YPoint{FT}, B}
const ZIntervalDomain{FT, B} = IntervalDomain{Geometry.ZPoint{FT}, B}

isperiodic(domain::IntervalDomain) = isnothing(domain.boundary_names)
boundary_names(domain::IntervalDomain) =
    isperiodic(domain) ? () : unique(domain.boundary_names)

"""
    IntervalDomain(coord⁻, coord⁺; periodic=true)
    IntervalDomain(coord⁻, coord⁺; boundary_names::Tuple{Symbol,Symbol})

Construct a `IntervalDomain`, the closed interval is given by `coord⁻`, `coord⁺` coordinate arguments.

Either a `periodic` or `boundary_names` keyword argument is required.
"""
function IntervalDomain(
    coord_min::Geometry.Abstract1DPoint,
    coord_max::Geometry.Abstract1DPoint;
    periodic = false,
    boundary_names::BCTagType = nothing,
)
    if !periodic && isnothing(boundary_names)
        throw(
            ArgumentError(
                "if `periodic=false` then an `boundary_names::Tuple{Symbol,Symbol}` keyword argument is required.",
            ),
        )
    end
    IntervalDomain(promote(coord_min, coord_max)..., boundary_names)
end
IntervalDomain(coords::ClosedInterval; kwargs...) =
    IntervalDomain(coords.left, coords.right; kwargs...)

"""
    z_max(domain::IntervalDomain)

The domain maximum along the z-direction.
"""
z_max(domain::IntervalDomain) = domain.coord_max.z

"""
    z_min(domain::IntervalDomain)

The domain minimum along the z-direction.
"""
z_min(domain::IntervalDomain) = domain.coord_min.z

function XIntervalDomain(;
    x_min::Real,
    x_max::Real,
    x_periodic::Bool = false,
    x_boundary_names = (:west, :east),
)
    x_min, x_max = promote(x_min, x_max)
    IntervalDomain(
        Geometry.XPoint(x_min),
        Geometry.XPoint(x_max);
        periodic = x_periodic,
        boundary_names = x_boundary_names,
    )
end
function YIntervalDomain(;
    y_min::Real,
    y_max::Real,
    y_periodic::Bool = false,
    y_boundary_names = (:south, :north),
)
    y_min, y_max = promote(y_min, y_max)
    IntervalDomain(
        Geometry.YPoint(y_min),
        Geometry.YPoint(y_max);
        periodic = y_periodic,
        boundary_names = y_boundary_names,
    )
end
function ZIntervalDomain(;
    z_min::Real,
    z_max::Real,
    z_periodic::Bool = false,
    z_boundary_names = (:bottom, :top),
)
    z_min, z_max = promote(z_min, z_max)
    IntervalDomain(
        Geometry.ZPoint(z_min),
        Geometry.ZPoint(z_max);
        periodic = z_periodic,
        boundary_names = z_boundary_names,
    )
end

coordinate_type(::IntervalDomain{CT}) where {CT} = CT
Base.eltype(domain::IntervalDomain) = coordinate_type(domain)

function print_interval(io::IO, domain::IntervalDomain{CT}) where {CT}
    print(
        io,
        fieldname(CT, 1),
        " ∈ [",
        Geometry.component(domain.coord_min, 1),
        ",",
        Geometry.component(domain.coord_max, 1),
        "] ",
    )
    if isperiodic(domain)
        print(io, "(periodic)")
    else
        print(io, domain.boundary_names)
    end
end
function Base.show(io::IO, domain::IntervalDomain)
    print(io, nameof(typeof(domain)), ": ")
    print_interval(io, domain)
end

struct RectangleDomain{I1 <: IntervalDomain, I2 <: IntervalDomain} <:
       AbstractDomain
    interval1::I1
    interval2::I2
end

const XYRectangleDomain = RectangleDomain{<:XIntervalDomain, <:YIntervalDomain}

Base.:*(interval1::IntervalDomain, interval2::IntervalDomain) =
    RectangleDomain(interval1, interval2)

boundary_names(domain::RectangleDomain) = unique(
    Symbol[
        boundary_names(domain.interval1)...,
        boundary_names(domain.interval2)...,
    ],
)::Vector{Symbol}

"""
    RectangleDomain(x1::ClosedInterval, x2::ClosedInterval;
        x1boundary::Tuple{Symbol,Symbol},
        x2boundary::Tuple{Symbol,Symbol},
        x1periodic = false,
        x2periodic = false,
    )

Construct a `RectangularDomain` in the horizontal.
If a given x1 or x2 boundary is not periodic, then `x1boundary` or `x2boundary` boundary name keyword arguments must be supplied.
"""
function RectangleDomain(
    x1::ClosedInterval{X1CT},
    x2::ClosedInterval{X2CT};
    x1periodic = false,
    x2periodic = false,
    x1boundary::BCTagType = nothing,
    x2boundary::BCTagType = nothing,
) where {X1CT <: Geometry.Abstract1DPoint, X2CT <: Geometry.Abstract1DPoint}
    interval1 =
        IntervalDomain(x1; periodic = x1periodic, boundary_names = x1boundary)
    interval2 =
        IntervalDomain(x2; periodic = x2periodic, boundary_names = x2boundary)
    return interval1 * interval2
end

function XYRectangleDomain(;
    x_min::Real,
    x_max::Real,
    x_periodic::Bool = false,
    x_boundary_names = (:west, :east),
    y_min::Real,
    y_max::Real,
    y_periodic::Bool = false,
    y_boundary_names = (:south, :north),
)
    x_domain = XIntervalDomain(; x_min, x_max, x_periodic, x_boundary_names)
    y_domain = YIntervalDomain(; y_min, y_max, y_periodic, y_boundary_names)
    RectangleDomain(x_domain, y_domain)
end

function Base.show(io::IO, domain::RectangleDomain)
    print(io, nameof(typeof(domain)), ": ")
    print_interval(io, domain.interval1)
    print(io, " × ")
    print_interval(io, domain.interval2)
end

coordinate_type(domain::RectangleDomain) = typeof(
    Geometry.product_coordinates(
        domain.interval1.coord_min,
        domain.interval2.coord_min,
    ),
)

"""
    SphereDomain(radius)

A domain representing the surface of a sphere with radius `radius`.
"""
struct SphereDomain{FT} <: AbstractDomain where {FT <: AbstractFloat}
    radius::FT
end
Base.show(io::IO, domain::SphereDomain) =
    print(io, nameof(typeof(domain)), ": radius = ", domain.radius)

boundary_names(::SphereDomain) = ()
coordinate_type(::SphereDomain{FT}) where {FT} = Geometry.Cartesian123Point{FT}

end # module