struct Line{PT}
    p::PT
    tangent::PT

    Line{PT}(p::PT, tangent::PT) where PT = new{PT}(p,tangent)
end

"""
    Line(p,t)

A line, as a callable object, starting at `p` with tangent `t`.

# Example
```julia-repl
julia> l = Grids.Line([1,1],[2,1])
Diffinitive.Grids.Line{StaticArraysCore.SVector{2, Int64}}([1, 1], [2, 1])

julia> l(0)
2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
 1
 1

julia> l(1)
2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
 3
 2
```

See also: [`LineSegment`](@ref).
"""
function Line(p, t)
    p = SVector{length(p)}(p)
    t = SVector{length(t)}(t)
    p, t = promote(p, t)

    return Line{typeof(p)}(p,t)
end

(c::Line)(s) = c.p + s*c.tangent

Grids.jacobian(l::Line, t) = l.tangent

struct LineSegment{PT}
    a::PT
    b::PT

    LineSegment{PT}(p::PT, tangent::PT) where PT = new{PT}(p,tangent)
end

"""
    LineSegment(a,b)

A line segment, as a callable object, from `a` to `b`.

# Example
```julia-repl
julia> l = Grids.LineSegment([1,1],[2,1])
Diffinitive.Grids.LineSegment{StaticArraysCore.SVector{2, Int64}}([1, 1], [2, 1])

julia> l(0)
2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
 1
 1

julia> l(0.5)
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 1.5
 1.0

julia> l(1)
2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
 2
 1
```

See also: [`Line`](@ref).
"""
function LineSegment(a, b)
    a = SVector{length(a)}(a)
    b = SVector{length(b)}(b)
    a, b = promote(a, b)

    return LineSegment{typeof(a)}(a,b)
end

(c::LineSegment)(s) = (1-s)*c.a + s*c.b

Grids.jacobian(c::LineSegment, s) = c.b - c.a

"""
    linesegments(ps...)

An array of line segments between the points `ps[1]`, `ps[2]`, and so on.

See also: [`polygon_edges`](@ref).
"""
function linesegments(ps...)
    return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1]
end


"""
    polygon_edges(ps...)

An array of line segments between the points `ps[1]`, `ps[2]`, and so on
including the segment between `ps[end]` and `ps[1]`.

See also: [`linesegments`](@ref).
"""
function polygon_edges(ps...)
    n = length(ps)
    return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(ps)]
end

struct Circle{PT,T}
    c::PT
    r::T

    Circle{PT,T}(c,r) where {PT,T} = new{PT,T}(c,r)
end

"""
    Circle(c,r)

A circle with center `c` and radius `r` paramatrized with the angle to the x-axis.

# Example
```julia-repl
julia> c = Grids.Circle([1,1], 2)
Diffinitive.Grids.Circle{StaticArraysCore.SVector{2, Int64}, Int64}([1, 1], 2)

julia> c(0)
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 3.0
 1.0

julia> c(π/2)
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 1.0000000000000002
 3.0
```
"""
function Circle(c,r)
    c = SVector{2}(c)
    return Circle{typeof(c), typeof(r)}(c,r)
end

function (C::Circle)(θ)
    (;c, r) = C
    c + r*@SVector[cos(θ), sin(θ)]
end

function Grids.jacobian(C::Circle, θ)
    (;r) = C
    r*@SVector[-sin(θ), cos(θ)]
end

"""
    TransfiniteInterpolationSurface(c₁, c₂, c₃, c₄)

A surface defined by the transfinite interpolation of the curves `c₁`, `c₂`, `c₃`, and  `c₄`.
"""
struct TransfiniteInterpolationSurface{T1,T2,T3,T4}
    c₁::T1
    c₂::T2
    c₃::T3
    c₄::T4
end

function (s::TransfiniteInterpolationSurface)(u,v)
    c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
    P₀₀ = c₁(0)
    P₁₀ = c₂(0)
    P₁₁ = c₃(0)
    P₀₁ = c₄(0)
    return (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - (
        (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁
    )
end

function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray)
    s(ξ̄...)
end

"""
    check_transfiniteinterpolation(s::TransfiniteInterpolationSurface)

Throw an error if the ends of the curves in the transfinite interpolation do not match.
"""
function check_transfiniteinterpolation(s::TransfiniteInterpolationSurface)
    if check_transfiniteinterpolation(Bool, s)
        return nothing
    else
        throw(ArgumentError("The end of each curve in the transfinite interpolation should be the same as the beginning of the next curve."))
    end
end

"""
    check_transfiniteinterpolation(::Type{Bool}, s::TransfiniteInterpolationSurface)

Return true if the ends of the curves in the transfinite interpolation match.
"""
function check_transfiniteinterpolation(::Type{Bool}, s::TransfiniteInterpolationSurface)
    if !isapprox(s.c₁(1), s.c₂(0))
        return false
    end

    if !isapprox(s.c₂(1), s.c₃(0))
        return false
    end

    if !isapprox(s.c₃(1), s.c₄(0))
        return false
    end

    if !isapprox(s.c₄(1), s.c₁(0))
        return false
    end

    return true
end

function Grids.jacobian(s::TransfiniteInterpolationSurface, ξ̄)
    u, v = ξ̄

    c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
    P₀₀ = c₁(0)
    P₁₀ = c₂(0)
    P₁₁ = c₃(0)
    P₀₁ = c₄(0)

    ∂x̄∂ξ₁ = (1-v)*jacobian(c₁,u) + c₂(v) - v*jacobian(c₃,1-u) -c₄(1-v) - (
        -(1-v)*P₀₀ + (1-v)*P₁₀ + v*P₁₁ - v*P₀₁
    )

    ∂x̄∂ξ₂ = -c₁(u) + u*jacobian(c₂,v) + c₃(1-u) - (1-u)*jacobian(c₄,1-v) - (
        -(1-u)*P₀₀ - u*P₁₀ + u*P₁₁ + (1-u)*P₀₁
    )

    return [∂x̄∂ξ₁ ∂x̄∂ξ₂]
end