changeset 2003:524a52f190d7 feature/sbp_operators/laplace_curvilinear

Merge feature/grids/geometry_functions
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Apr 2025 09:03:05 +0200
parents 496b8d3903c6 (current diff) 4300c59bbeff (diff)
children
files Manifest.toml src/Grids/Grids.jl
diffstat 16 files changed, 617 insertions(+), 74 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Tue Feb 11 09:03:04 2025 +0100
+++ b/Manifest.toml	Tue Apr 29 09:03:05 2025 +0200
@@ -60,9 +60,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01"
+git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.11"
+version = "1.9.12"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/benchmark/Manifest.toml	Tue Feb 11 09:03:04 2025 +0100
+++ b/benchmark/Manifest.toml	Tue Apr 29 09:03:05 2025 +0200
@@ -61,7 +61,7 @@
 deps = ["LinearAlgebra", "StaticArrays", "TOML"]
 path = ".."
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-version = "0.1.3"
+version = "0.1.4"
 
     [deps.Diffinitive.extensions]
     DiffinitiveMakieExt = "Makie"
@@ -242,9 +242,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01"
+git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.11"
+version = "1.9.12"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/docs/Manifest.toml	Tue Feb 11 09:03:04 2025 +0100
+++ b/docs/Manifest.toml	Tue Apr 29 09:03:05 2025 +0200
@@ -28,9 +28,9 @@
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759"
+git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.6"
+version = "0.7.8"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -46,7 +46,7 @@
 deps = ["LinearAlgebra", "StaticArrays", "TOML"]
 path = ".."
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-version = "0.1.3"
+version = "0.1.4"
 
     [deps.Diffinitive.extensions]
     DiffinitiveMakieExt = "Makie"
@@ -67,9 +67,9 @@
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
-git-tree-sha1 = "d0ea2c044963ed6f37703cead7e29f70cba13d7e"
+git-tree-sha1 = "182a9a3fe886587ba230a417f1651a4cbc2b92d4"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.8.0"
+version = "1.8.1"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -78,9 +78,9 @@
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e51db81749b0777b2147fbe7b783ee79045b8e99"
+git-tree-sha1 = "d55dffd9ae73ff72f1c0482454dcf2ec6c6c4a63"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.6.4+3"
+version = "2.6.5+0"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
@@ -205,9 +205,9 @@
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10"
+git-tree-sha1 = "a9697f1d06cc3eb3fb3ad49cc67f2cfabaac31ea"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.0.15+3"
+version = "3.0.16+0"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -276,9 +276,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01"
+git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.11"
+version = "1.9.12"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/docs/make.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/docs/make.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -30,7 +30,7 @@
     "operator_file_format.md",
     "grids_and_grid_functions.md",
     "matrix_and_tensor_representations.md",
-    "manifolds_charts_atlases.md"
+    "manifolds_charts_atlases.md",
     "Submodules" => [
         "submodules/grids.md",
         "submodules/lazy_tensors.md",
--- a/docs/src/manifolds_charts_atlases.md	Tue Feb 11 09:03:04 2025 +0100
+++ b/docs/src/manifolds_charts_atlases.md	Tue Apr 29 09:03:05 2025 +0200
@@ -1,18 +1,25 @@
 # Manifolds, Charts, and Atlases
 
 To construct grids on more complicated geometries we use manifolds described
-by one or more charts. The charts describe a mapping from some parameter space
-to the geometry that we are interested in. If there are more than one chart
-for a given geometry this collection of charts and their connection is
-described by and atlas.
+by one or more charts. The charts describe a mapping from a logical parameter
+space to the geometry that we are interested in. If there are more than one
+chart for a given geometry this collection of charts and how they are
+connected is described by an atlas.
+
+We consider a mapping from the logical coordidinates ``\xi \in \Xi`` to the
+physical coordinates ``x \in \Omega``. A `Chart` describes the mapping by a
+`ParameterSpace` respresenting ``\Xi`` and some mapping object that takes
+arguments ``\xi \in \Xi`` and returns coordinates ``x\in\Omega``. The mapping
+object can either be a function or some other callable object.
 
 For the construction of differential and difference operators on a manifold
-with a chart the library needs to know the Jacobian of the mapping as a
-function of coordinates in the logical parameter space. Internally,
-Diffinitive.jl uses a local Jacobian function, `Grids.jacobian(f, ξ)`. For
-geometry objects provided by the library this function should have fast and
-efficient implementations. If you are creating your own mapping functions you
-can implement `Grids.jacobian` for your function or type, for example
+with a chart the library needs to know the Jacobian,
+``\frac{\partial x}{\partial \xi}``, of the mapping as a function of
+coordinates in the logical parameter space. Internally, Diffinitive.jl uses a
+local Jacobian function, `Grids.jacobian(f, ξ)`. For geometry objects provided
+by the library this function should have fast and efficient implementations.
+If you are creating your own mapping functions you must implement
+`Grids.jacobian` for your function or type, for example
 
 ```julia
 f(x) = 2x
@@ -32,5 +39,3 @@
 using ForwardDiff
 Grids.jacobian(f,x) = ForwardDiff.jacobian(f,x)
 ```
-
-<!-- What more needs to be said here? --/>
--- a/ext/DiffinitivePlotsExt.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/ext/DiffinitivePlotsExt.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -57,16 +57,4 @@
     return x, y
 end
 
-# get_axis_limits(plt, :x)
-
-
-# ReicpesPipline/src/user_recipe.jl
-# @recipe function f(f::FuncOrFuncs{F}) where {F<:Function}
-
-# @recipe function f(f::Function, xmin::Number, xmax::Number)
-
-# _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
-
 end
-
-
--- a/src/Grids/Grids.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/src/Grids/Grids.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -84,4 +84,17 @@
 include("mapped_grid.jl")
 include("geometry.jl")
 
+function __init__()
+    if !isdefined(Base.Experimental, :register_error_hint)
+        return
+    end
+
+    Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
+        if exc.f == Grids.jacobian
+            print(io, "\nThis possibly means that a function used to define a coordinate mapping is missing a method for `Grids.jacobian`.\n")
+            print(io, "Provide one by for exmple implementing `Grids.jacobian(::$(typeof(exc.args[1])), x) = ...` or `Grids.jacobian(f, x) = ForwardDiff.jacobian(f,x)`")
+        end
+    end
+end
+
 end # module
--- a/src/Grids/geometry.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/src/Grids/geometry.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -1,32 +1,148 @@
 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{T,PT}
+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)(θ)
@@ -34,6 +150,16 @@
     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
@@ -56,4 +182,60 @@
     s(ξ̄...)
 end
 
-# TODO: Implement jacobian() for the different mapping helpers
+"""
+    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(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
--- a/src/Grids/manifolds.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/src/Grids/manifolds.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -9,9 +9,15 @@
 end
 
 Base.ndims(::Chart{D}) where D = D
-(c::Chart)(ξ) = c.mapping(ξ)
 parameterspace(c::Chart) = c.parameterspace
 
+function (c::Chart)(ξ)
+    if ξ ∉ parameterspace(c)
+        throw(DomainError(ξ, "chart was called logical coordinates outside the parameterspace. If this was inteded, use the `mapping` field from the Chart struct instead."))
+    end
+    return c.mapping(ξ)
+end
+
 """
     jacobian(c::Chart, ξ)
 
@@ -30,8 +36,12 @@
 ```
 which will both allow calling `jacobian(c,ξ)`.
 """
-jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ)
-# TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian?
+function jacobian(c::Chart, ξ)
+    if ξ ∉ parameterspace(c)
+        throw(DomainError(ξ, "jacobian was called with logical coordinates outside the parameterspace of the chart. If this was inteded, use the `mapping` field from the Chart struct instead."))
+    end
+    return jacobian(c.mapping, ξ)
+end
 
 boundary_identifiers(c::Chart) = boundary_identifiers(parameterspace(c))
 
@@ -54,7 +64,7 @@
 """
     connections(::Atlas)
 
-Collection of pairs of multiblock boundary identifiers.
+Collection of 2-tuples of multiblock boundary identifiers.
 """
 function connections end
 
--- a/src/Grids/parameter_space.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/src/Grids/parameter_space.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -18,6 +18,13 @@
 abstract type ParameterSpace{D} end
 Base.ndims(::ParameterSpace{D}) where D = D
 
+@doc """
+    in(x, S::ParameterSpace)
+    ∈(x, S::ParameterSpace)
+
+Test if the point `x` is in the parameter space `S`.
+""" Base.in(x,::ParameterSpace)
+
 """
     Interval{T} <: ParameterSpace{1}
 
@@ -47,6 +54,8 @@
 
 boundary_identifiers(::Interval) = (LowerBoundary(), UpperBoundary())
 
+Base.in(x, i::Interval) = i.a <= x <= i.b
+
 """
     unitinterval(T=Float64)
 
@@ -102,6 +111,11 @@
     end
 end
 
+function Base.in(x, box::HyperBox)
+    return all(eachindex(x)) do i
+        box.a[i] <= x[i] <= box.b[i]
+    end
+end
 
 """
     unitsquare(T=Float64)
@@ -150,6 +164,20 @@
     return Simplex(Tuple(convert(T,v) for v ∈ verticies))
 end
 
+function Base.in(x, s::Simplex)
+    v₁ = s.verticies[1]
+    V = map(s.verticies) do v
+        v - v₁
+    end
+
+    A = hcat(V[2:end]...) # Matrix with edge vectors as columns
+    λ = A \ (x - v₁)
+
+    λ_full = (1 - sum(λ), λ...) # Full barycentric coordinates
+
+    return all(λᵢ -> zero(λᵢ) ≤ λᵢ ≤ one(λᵢ), λ_full)
+end
+
 """
     verticies(s::Simplex)
 
--- a/test/Grids/equidistant_grid_test.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/test/Grids/equidistant_grid_test.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -164,7 +164,6 @@
         @test equidistant_grid(HyperBox((0,),(2,)),4) == equidistant_grid(@SVector[0], @SVector[2], 4)
     end
 
-
     @testset "equidistant_grid(::Chart)" begin
         c = Chart(unitsquare()) do (ξ,η)
             @SVector[2ξ, 3η]
--- a/test/Grids/geometry_test.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/test/Grids/geometry_test.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -1,23 +1,274 @@
+using Diffinitive.Grids
+using Diffinitive.Grids: Line, LineSegment, linesegments, polygon_edges, Circle, TransfiniteInterpolationSurface, check_transfiniteinterpolation
+using StaticArrays
+
 @testset "Line" begin
-    @test_broken false
+    @testset "Constructors" begin
+        @test Line([1,2],[2,3]) isa Line{SVector{2,Int}}
+        @test Line((1,2),(2,3)) isa Line{SVector{2,Int}}
+        @test Line(@SVector[1,2],[2,3]) isa Line{SVector{2,Int}}
+        @test Line(@SVector[1,2],@SVector[2,3]) isa Line{SVector{2,Int}}
+
+        @test Line([1,2],[2.,3]) isa Line{SVector{2,Float64}}
+        @test Line(@SVector[1,2.],@SVector[2,3]) isa Line{SVector{2,Float64}}
+        @test Line((1,2.),(2,3)) isa Line{SVector{2,Float64}}
+    end
+
+    @testset "Evaluation" begin
+        l = Line([1,2],[2,3])
+
+        @test l(0) == [1,2]
+        @test l(1) == [1,2] + [2,3]
+        @test l(1/2) == [1,2] + [2,3]/2
+    end
+
+    @testset "Grids.jacobian" begin
+        l = Line([1,2],[2,3])
+
+        @test Grids.jacobian(l,0) == [2,3]
+        @test Grids.jacobian(l,1) == [2,3]
+        @test Grids.jacobian(l,1/2) == [2,3]
+    end
 end
 
 @testset "LineSegment" begin
-    @test_broken false
+    @testset "Constructors" begin
+        @test LineSegment([1,2],[2,3]) isa LineSegment{SVector{2,Int}}
+        @test LineSegment((1,2),(2,3)) isa LineSegment{SVector{2,Int}}
+        @test LineSegment(@SVector[1,2],[2,3]) isa LineSegment{SVector{2,Int}}
+        @test LineSegment(@SVector[1,2],@SVector[2,3]) isa LineSegment{SVector{2,Int}}
+
+        @test LineSegment([1,2],[2.,3]) isa LineSegment{SVector{2,Float64}}
+        @test LineSegment(@SVector[1,2.],@SVector[2,3]) isa LineSegment{SVector{2,Float64}}
+        @test LineSegment((1,2.),(2,3)) isa LineSegment{SVector{2,Float64}}
+    end
+
+    @testset "Evaluation" begin
+        l = LineSegment([1,2],[2,3])
+
+        @test l(0) == [1,2]
+        @test l(1) == [2,3]
+        @test l(1/2) == [1,2]/2 + [2,3]/2
+    end
+
+    @testset "Grids.jacobian" begin
+        a, b = [1,2], [2,3]
+        l = LineSegment(a,b)
+        d = b-a
+
+        @test Grids.jacobian(l,0) == d
+        @test Grids.jacobian(l,1) == d
+        @test Grids.jacobian(l,1/2) == d
+    end
 end
 
 @testset "linesegments" begin
-    @test_broken false
+    a,b,c,d = [1,1],[2,2],[3,3],[4,4]
+    @test linesegments(a,b) == [
+        LineSegment(a,b),
+    ]
+
+    @test linesegments(a,b,c) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+    ]
+
+    @test linesegments(a,b,c,d) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+        LineSegment(c,d),
+    ]
 end
 
 @testset "polygon_edges" begin
-    @test_broken false
+    a,b,c,d = [1,1],[2,2],[3,3],[4,4]
+    @test polygon_edges(a,b) == [
+        LineSegment(a,b),
+        LineSegment(b,a),
+    ]
+
+    @test polygon_edges(a,b,c) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+        LineSegment(c,a),
+    ]
+
+    @test polygon_edges(a,b,c,d) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+        LineSegment(c,d),
+        LineSegment(d,a),
+    ]
 end
 
 @testset "Circle" begin
-    @test_broken false
+    @testset "Constructors" begin
+        @test Circle([1,2], 1) isa Circle{SVector{2,Int},Int}
+        @test Circle([1,2], 1.) isa Circle{SVector{2,Int},Float64}
+        @test Circle([1,2.], 1.) isa Circle{SVector{2,Float64},Float64}
+        @test Circle([1,2.], 1) isa Circle{SVector{2,Float64},Int}
+        @test Circle((1,2.), 1.) isa Circle{SVector{2,Float64},Float64}
+        @test Circle((1,2), 1.) isa Circle{SVector{2,Int},Float64}
+        @test Circle((1.,2), 1) isa Circle{SVector{2,Float64},Int}
+        @test Circle((1,2), 1) isa Circle{SVector{2,Int},Int}
+        @test Circle(@SVector[1,2], 1.) isa Circle{SVector{2,Int},Float64}
+        @test Circle(@SVector[1,2.], 1.) isa Circle{SVector{2,Float64},Float64}
+    end
+
+    @testset "Evaluation" begin
+        c = Circle([0,0], 1)
+        @test c(0) ≈ [1,0]
+        @test c(π/2) ≈ [0,1]
+        @test c(π) ≈ [-1,0]
+        @test c(3π/2) ≈ [0,-1]
+        @test c(π/4) ≈ [1/√(2),1/√(2)]
+
+        c = Circle([0,0], 2)
+        @test c(0) ≈ [2,0]
+        @test c(π/2) ≈ [0,2]
+        @test c(π) ≈ [-2,0]
+        @test c(3π/2) ≈ [0,-2]
+        @test c(π/4) ≈ [√(2),√(2)]
+    end
+
+    @testset "Grids.jacobian" begin
+        c = Circle([0,0], 1)
+        @test Grids.jacobian(c, 0) ≈ [0,1]
+        @test Grids.jacobian(c, π/2) ≈ [-1,0]
+        @test Grids.jacobian(c, π) ≈ [0,-1]
+        @test Grids.jacobian(c, 3π/2) ≈ [1,0]
+        @test Grids.jacobian(c, π/4) ≈ [-1/√(2),1/√(2)]
+
+        c = Circle([0,0], 2)
+        @test Grids.jacobian(c, 0) ≈ [0,2]
+        @test Grids.jacobian(c, π/2) ≈ [-2,0]
+        @test Grids.jacobian(c, π) ≈ [0,-2]
+        @test Grids.jacobian(c, 3π/2) ≈ [2,0]
+        @test Grids.jacobian(c, π/4) ≈ [-√(2),√(2)]
+
+        c = Circle([-1,1], 1)
+        @test Grids.jacobian(c, 0) ≈ [0,1]
+        @test Grids.jacobian(c, π/2) ≈ [-1,0]
+        @test Grids.jacobian(c, π) ≈ [0,-1]
+        @test Grids.jacobian(c, 3π/2) ≈ [1,0]
+        @test Grids.jacobian(c, π/4) ≈ [-1/√(2),1/√(2)]
+
+        c = Circle([-1,1], 2)
+        @test Grids.jacobian(c, 0) ≈ [0,2]
+        @test Grids.jacobian(c, π/2) ≈ [-2,0]
+        @test Grids.jacobian(c, π) ≈ [0,-2]
+        @test Grids.jacobian(c, 3π/2) ≈ [2,0]
+        @test Grids.jacobian(c, π/4) ≈ [-√(2),√(2)]
+    end
 end
 
 @testset "TransfiniteInterpolationSurface" begin
-    @test_broken false
+    @testset "Constructors" begin
+        @test TransfiniteInterpolationSurface(t->[1,2], t->[2,1], t->[0,0], t->[1,1]) isa TransfiniteInterpolationSurface
+
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        @test TransfiniteInterpolationSurface(cs...) isa TransfiniteInterpolationSurface
+    end
+
+    @testset "Evaluation" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test ti(0,0) == [0,0]
+        @test ti([0,0]) == [0,0]
+        @test ti(1,0) == [1,0]
+        @test ti([1,0]) == [1,0]
+        @test ti(1,1) == [1,1]
+        @test ti([1,1]) == [1,1]
+        @test ti(0,1) == [0,1]
+        @test ti([0,1]) == [0,1]
+
+        @test ti(1/2, 0) == [1/2, 0]
+        @test ti(1/2, 1) == [1/2, 1]
+        @test ti(0,1/2) == [0,1/2]
+        @test ti(1,1/2) == [1,1/2]
+
+
+        a, b, c, d = [1,0],[2,1/4],[2.5,1],[-1/3,1]
+        cs = polygon_edges(a,b,c,d)
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test ti(0,0) == a
+        @test ti(1,0) == b
+        @test ti(1,1) == c
+        @test ti(0,1) == d
+
+        @test ti(1/2, 0) == (a+b)/2
+        @test ti(1/2, 1) == (c+d)/2
+        @test ti(0, 1/2) == (a+d)/2
+        @test ti(1, 1/2) == (b+c)/2
+
+        # TODO: Some test with curved edges?
+    end
+
+    @testset "check_transfiniteinterpolation" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test check_transfiniteinterpolation(ti) == nothing
+        @test check_transfiniteinterpolation(Bool, ti) == true
+
+        bad_sides = [
+            LineSegment([0,0],[1/2,0]),
+            LineSegment([1,0],[1,1/2]),
+            LineSegment([1,1],[0,1/2]),
+            LineSegment([0,1],[0,1/2]),
+        ]
+
+        s1 = TransfiniteInterpolationSurface(bad_sides[1],cs[2],cs[3],cs[4])
+        s2 = TransfiniteInterpolationSurface(cs[1],bad_sides[2],cs[3],cs[4])
+        s3 = TransfiniteInterpolationSurface(cs[1],cs[2],bad_sides[3],cs[4])
+        s4 = TransfiniteInterpolationSurface(cs[1],cs[2],cs[3],bad_sides[4])
+
+        @test check_transfiniteinterpolation(Bool, s1) == false
+        @test check_transfiniteinterpolation(Bool, s2) == false
+        @test check_transfiniteinterpolation(Bool, s3) == false
+        @test check_transfiniteinterpolation(Bool, s4) == false
+
+        @test_throws ArgumentError check_transfiniteinterpolation(s1)
+        @test_throws ArgumentError check_transfiniteinterpolation(s2)
+        @test_throws ArgumentError check_transfiniteinterpolation(s3)
+        @test_throws ArgumentError check_transfiniteinterpolation(s4)
+    end
+
+    @testset "Grids.jacobian" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test Grids.jacobian(ti, [0,0]) isa SMatrix
+
+        @test Grids.jacobian(ti, [0,0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,1]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [0,1]) == [1 0; 0 1]
+
+        @test Grids.jacobian(ti, [1/2, 0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1/2, 1]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [0,1/2]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,1/2]) == [1 0; 0 1]
+
+
+        a, b, c, d = [1,0],[2,1/4],[2.5,1],[-1/3,1]
+        cs = polygon_edges(a,b,c,d)
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test Grids.jacobian(ti, [0,0]) == [b-a d-a]
+        @test Grids.jacobian(ti, [1,0]) == [b-a c-b]
+        @test Grids.jacobian(ti, [1,1]) == [c-d c-b]
+        @test Grids.jacobian(ti, [0,1]) == [c-d d-a]
+
+
+        mid(x,y) = (x+y)/2
+        @test Grids.jacobian(ti, [1/2, 0]) ≈ [b-a mid(c,d)-mid(a,b)]
+        @test Grids.jacobian(ti, [1/2, 1]) ≈ [c-d mid(c,d)-mid(a,b)]
+        @test Grids.jacobian(ti, [0, 1/2]) ≈ [mid(b,c)-mid(a,d) d-a]
+        @test Grids.jacobian(ti, [1, 1/2]) ≈ [mid(b,c)-mid(a,d) c-b]
+
+        # TODO: Some test with curved edges?
+    end
 end
--- a/test/Grids/manifolds_test.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/test/Grids/manifolds_test.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -15,14 +15,23 @@
 
 @testset "Chart" begin
     X(ξ) = 2ξ
-    Grids.jacobian(::typeof(X), ξ) = @SVector[2,2]
+    Grids.jacobian(::typeof(X), ξ) = @SMatrix[2 0; 0 2]
     c = Chart(X, unitsquare())
     @test c isa Chart{2}
-    @test c([3,2]) == [6,4]
     @test parameterspace(c) == unitsquare()
     @test ndims(c) == 2
 
-    @test jacobian(c, [3,2]) == [2,2]
+    @testset "Calling" begin
+        c = Chart(X, unitsquare())
+        @test c([.3,.2]) == [.6,.4]
+        @test_throws DomainError c([3,2])
+    end
+
+    @testset "jacobian(::Chart, ⋅)" begin
+        c = Chart(X, unitsquare())
+        @test_throws DomainError jacobian(c, [3,2])
+        @test jacobian(c, [.3,.2]) == [2 0; 0 2]
+    end
 
     @test Set(boundary_identifiers(Chart(X,unitsquare()))) == Set([east(),west(),south(),north()])
 end
--- a/test/Grids/parameter_space_test.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/test/Grids/parameter_space_test.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -1,6 +1,7 @@
 using Test
 
 using Diffinitive.Grids
+using StaticArrays
 
 @testset "ParameterSpace" begin
     @test ndims(HyperBox([1,1], [2,2])) == 2
@@ -22,6 +23,13 @@
     @test limits(unitinterval(Int)) == (0,1)
 
     @test boundary_identifiers(unitinterval()) == (LowerBoundary(), UpperBoundary())
+
+    @test 0 ∈ Interval(0,1)
+    @test 0. ∈ Interval(0,1)
+    @test 1. ∈ Interval(0,1)
+    @test 2 ∉ Interval(0,1)
+    @test -1 ∉ Interval(0,1)
+    @test -1. ∉ Interval(0,1)
 end
 
 @testset "HyperBox" begin
@@ -59,6 +67,13 @@
         CartesianBoundary{3,LowerBoundary}(),
         CartesianBoundary{3,UpperBoundary}(),
     ]
+
+    @test @SVector[1.5, 3.5] ∈ HyperBox([1,2], [3,4])
+    @test @SVector[1, 2] ∈ HyperBox([1,2], [3,4])
+    @test @SVector[3, 4] ∈ HyperBox([1,2], [3,4])
+
+    @test @SVector[0.5, 3.5] ∉ HyperBox([1,2], [3,4])
+    @test @SVector[1.5, 4.5] ∉ HyperBox([1,2], [3,4])
 end
 
 @testset "Simplex" begin
@@ -77,4 +92,32 @@
     @test verticies(unittetrahedron()) == ([0,0,0], [1,0,0], [0,1,0],[0,0,1])
 
     @test unitsimplex(4) isa Simplex{Float64,4}
+
+    @testset "Base.in" begin
+        @testset "2D" begin
+            T₂ = Simplex([0.0, 0.0], [1.0, 0.0], [0.0, 1.0])
+            @test [0.1, 0.1] ∈ T₂
+            @test [0.3, 0.3] ∈ T₂
+            @test [1.0, 0.0] ∈ T₂
+            @test [0.0, 0.0] ∈ T₂
+            @test [0.0, 1.0] ∈ T₂
+            @test [0.5, 0.5] ∈ T₂
+
+            @test [0.6, 0.6] ∉ T₂
+            @test [-0.1, 0.1] ∉ T₂
+        end
+
+        @testset "3D" begin
+            T₃ = Simplex([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0])
+            @test [0.1, 0.1, 0.1] ∈ T₃
+            @test [0.0, 0.0, 0.0] ∈ T₃
+            @test [1.0, 0.0, 0.0] ∈ T₃
+            @test [0.25, 0.25, 0.25] ∈ T₃
+            @test [0.5, 0.5, 0.0] ∈ T₃
+            @test [0.3, 0.3, 0.3] ∈ T₃
+
+            @test [0.5, 0.5, 1.0] ∉ T₃
+            @test [0.3, 0.3, 0.5] ∉ T₃
+        end
+    end
 end
--- a/test/Manifest.toml	Tue Feb 11 09:03:04 2025 +0100
+++ b/test/Manifest.toml	Tue Apr 29 09:03:05 2025 +0200
@@ -6,9 +6,9 @@
 
 [[deps.Aqua]]
 deps = ["Compat", "Pkg", "Test"]
-git-tree-sha1 = "3f96fac9ed31d16aba9f2f46fb5cbfc98db6b57c"
+git-tree-sha1 = "500941611ff835a025f484f55836f6feea61720a"
 uuid = "4c88cf16-eb10-579e-8560-4a9242c79595"
-version = "0.8.10"
+version = "0.8.11"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -101,14 +101,14 @@
 version = "0.2.4"
 
 [[deps.JET]]
-deps = ["CodeTracking", "InteractiveUtils", "JuliaInterpreter", "LoweredCodeUtils", "MacroTools", "Pkg", "PrecompileTools", "Preferences", "Test"]
-git-tree-sha1 = "24bdbf3ef611b69d1f5ef9331e69b6007152989e"
+deps = ["CodeTracking", "InteractiveUtils", "JuliaInterpreter", "JuliaSyntax", "LoweredCodeUtils", "MacroTools", "Pkg", "PrecompileTools", "Preferences", "Test"]
+git-tree-sha1 = "a453c9b3320dd73f5b05e8882446b6051cb254c4"
 uuid = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
-version = "0.9.14"
+version = "0.9.18"
 
     [deps.JET.extensions]
     JETCthulhuExt = "Cthulhu"
-    ReviseExt = "Revise"
+    JETReviseExt = "Revise"
 
     [deps.JET.weakdeps]
     Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
@@ -128,9 +128,14 @@
 
 [[deps.JuliaInterpreter]]
 deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
-git-tree-sha1 = "a729439c18f7112cbbd9fcdc1771ecc7f071df6a"
+git-tree-sha1 = "4bf4b400a8234cff0f177da4a160a90296159ce9"
 uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
-version = "0.9.39"
+version = "0.9.41"
+
+[[deps.JuliaSyntax]]
+git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a"
+uuid = "70703baa-626e-46a2-a12c-08ffd08c73b4"
+version = "0.4.10"
 
 [[deps.LRUCache]]
 git-tree-sha1 = "b3cc6698599b10e652832c2f23db3cab99d51b59"
@@ -226,9 +231,9 @@
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
-git-tree-sha1 = "fe891aea7ccd23897520db7f16931212454e277e"
+git-tree-sha1 = "cc0a5deefdb12ab3a096f00a6d42133af4560d71"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "1.1.1"
+version = "1.1.2"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -351,9 +356,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01"
+git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.11"
+version = "1.9.12"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/update_manifest.jl	Tue Feb 11 09:03:04 2025 +0100
+++ b/update_manifest.jl	Tue Apr 29 09:03:05 2025 +0200
@@ -1,12 +1,22 @@
 using Pkg
 
+function update_and_log(root, d)
+    printstyled("Updating '$d'"; color=:cyan, bold=true)
+    println()
+    dp = joinpath(root, d)
+
+    update_directory(dp)
+    println()
+end
+
 function update_directory(d)
     Pkg.activate(d)
     Pkg.update()
-    println()
 end
 
-update_directory(".")
-update_directory("benchmark")
-update_directory("docs")
-update_directory("test")
+root = @__DIR__
+update_and_log(root, ".")
+update_and_log(root, "benchmark")
+update_and_log(root, "docs")
+update_and_log(root, "test")
+