--- a/Manifest.toml	Tue Apr 29 08:48:58 2025 +0200
+++ b/Manifest.toml	Tue Apr 29 08:59:29 2025 +0200
@@ -2,7 +2,7 @@
 
 julia_version = "1.11.3"
 manifest_format = "2.0"
-project_hash = "383b10bc09a08d519764abf167e40dae29b16b2e"
+project_hash = "51374910f6b7831392791bcca0079860f1023f47"
 
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

--- a/Project.toml	Tue Apr 29 08:48:58 2025 +0200
+++ b/Project.toml	Tue Apr 29 08:59:29 2025 +0200
@@ -9,6 +9,7 @@
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 
 [weakdeps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -16,6 +17,7 @@
 
 [extensions]
 DiffinitiveMakieExt = "Makie"
+DiffinitivePlotsExt = "Plots"
 DiffinitiveSparseArrayKitExt = ["SparseArrayKit", "Tokens"]
 DiffinitiveSparseArraysExt = ["SparseArrays", "Tokens"]
 
@@ -24,6 +26,7 @@
 StaticArrays = "1.0"
 TOML = "1.0"
 Makie = "0.21"
+Plots = "1.40"
 SparseArrayKit = "0.4"
 Tokens = "0.1.1"
 SparseArrays = "1.10"

--- a/docs/make.jl	Tue Apr 29 08:48:58 2025 +0200
+++ b/docs/make.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -30,6 +30,7 @@
     "operator_file_format.md",
     "grids_and_grid_functions.md",
     "matrix_and_tensor_representations.md",
+    "manifolds_charts_atlases.md",
     "Submodules" => [
         "submodules/grids.md",
         "submodules/lazy_tensors.md",

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/manifolds_charts_atlases.md	Tue Apr 29 08:59:29 2025 +0200
@@ -0,0 +1,41 @@
+# 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 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,
+``\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
+Grids.jacobian(::typeof(f), x) = fill(2, length(x))
+```
+
+```julia
+struct F end
+(::F)(x) = 2x
+Grids.jacobian(::F, x) = fill(2,length(x))
+```
+
+You can also provide a fallback function using one of the many automatic
+differentiation packages, for example
+
+```julia
+using ForwardDiff
+Grids.jacobian(f,x) = ForwardDiff.jacobian(f,x)
+```

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/DiffinitivePlotsExt.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -0,0 +1,60 @@
+module DiffinitivePlotsExt
+
+using Diffinitive.Grids
+using Plots
+
+@recipe f(::Type{<:Grid}, g::Grid) = map(Tuple,g)[:]
+
+@recipe function f(c::Chart{2,<:Rectangle}, n=5, m=n; draw_border=true, bordercolor=1)
+    Ξ = parameterspace(c)
+    ξs = range(limits(Ξ,1)..., n)
+    ηs = range(limits(Ξ,2)..., m)
+
+    label := false
+    seriescolor --> 2
+    for ξ ∈ ξs
+        @series adapted_curve_grid(η->c((ξ,η)),limits(Ξ,1))
+    end
+
+    for η ∈ ηs
+        @series adapted_curve_grid(ξ->c((ξ,η)),limits(Ξ,2))
+    end
+
+    if ~draw_border
+        return
+    end
+
+    for ξ ∈ limits(Ξ,1)
+        @series begin
+            linewidth --> 3
+            seriescolor := bordercolor
+            adapted_curve_grid(η->c((ξ,η)),limits(Ξ,1))
+        end
+    end
+
+    for η ∈ limits(Ξ,2)
+        @series begin
+            linewidth --> 3
+            seriescolor := bordercolor
+            adapted_curve_grid(ξ->c((ξ,η)),limits(Ξ,2))
+        end
+    end
+end
+
+function adapted_curve_grid(g, minmax)
+    t1, _ = PlotUtils.adapted_grid(t->g(t)[1], minmax)
+    t2, _ = PlotUtils.adapted_grid(t->g(t)[2], minmax)
+
+    ts = sort(vcat(t1,t2))
+
+    x = map(ts) do t
+        g(t)[1]
+    end
+    y = map(ts) do t
+        g(t)[2]
+    end
+
+    return x, y
+end
+
+end

--- a/src/Grids/Grids.jl	Tue Apr 29 08:48:58 2025 +0200
+++ b/src/Grids/Grids.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -24,6 +24,16 @@
 export unittetrahedron
 export unitsimplex
 
+export Chart
+
+export Atlas
+export charts
+export connections
+export CartesianAtlas
+export UnstructuredAtlas
+
+export parameterspace
+
 # Grid
 export Grid
 export coordinate_size
@@ -66,10 +76,24 @@
 
 include("parameter_space.jl")
 include("grid.jl")
+include("multiblockgrids.jl")
+include("manifolds.jl")
 include("tensor_grid.jl")
 include("equidistant_grid.jl")
 include("zero_dim_grid.jl")
 include("mapped_grid.jl")
-include("multiblockgrids.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/equidistant_grid.jl	Tue Apr 29 08:48:58 2025 +0200
+++ b/src/Grids/equidistant_grid.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -154,6 +154,10 @@
 equidistant_grid(d::Interval, size::Int) = equidistant_grid(limits(d)..., size)
 equidistant_grid(hb::HyperBox, dims::Vararg{Int}) = equidistant_grid(limits(hb)..., dims...)
 
+function equidistant_grid(c::Chart, dims::Vararg{Int})
+    mapped_grid(c, ξ->jacobian(c,ξ), parameterspace(c), dims...)
+end
+
 
 CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/manifolds.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -0,0 +1,162 @@
+"""
+    Chart{D}
+
+A parametrized description of a manifold or part of a manifold.
+"""
+struct Chart{D, PST<:ParameterSpace{D}, MT}
+    mapping::MT
+    parameterspace::PST
+end
+
+Base.ndims(::Chart{D}) where D = D
+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, ξ)
+
+The jacobian of the mapping evaluated at `ξ`. This defers to the
+implementation of `jacobian` for the mapping itself. If no implementation is
+available one can easily be specified for either the mapping function or the
+chart itself.
+```julia
+c = Chart(f, ps)
+jacobian(f::typeof(f), ξ) = f′(ξ)
+```
+or
+```julia
+c = Chart(f, ps)
+jacobian(c::typeof(c),ξ) = f′(ξ)
+```
+which will both allow calling `jacobian(c,ξ)`.
+"""
+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))
+
+
+"""
+    Atlas
+
+A collection of charts and their connections.
+Should implement methods for `charts` and `connections`.
+"""
+abstract type Atlas end
+
+"""
+    charts(::Atlas)
+
+The colloction of charts in the atlas.
+"""
+function charts end
+
+"""
+    connections(::Atlas)
+
+Collection of 2-tuples of multiblock boundary identifiers.
+"""
+function connections end
+
+
+"""
+    CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas
+
+An atlas where the charts are arranged and connected like an array.
+"""
+struct CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas
+    charts::AT
+end
+
+charts(a::CartesianAtlas) = a.charts
+
+function connections(a::CartesianAtlas)
+    c = Tuple{MultiBlockBoundary, MultiBlockBoundary}[]
+
+    for d ∈ 1:ndims(charts(a))
+        Is = eachslice(CartesianIndices(charts(a)); dims=d)
+        for i ∈ 1:length(Is)-1 # For each interface between slices
+            for jk ∈ eachindex(Is[i]) # For each block in slice
+                Iᵢⱼₖ = Tuple(Is[i][jk])
+                Iᵢ₊₁ⱼₖ = Tuple(Is[i+1][jk])
+                push!(c,
+                    (
+                        MultiBlockBoundary{Iᵢⱼₖ,   CartesianBoundary{d,UpperBoundary}}(),
+                        MultiBlockBoundary{Iᵢ₊₁ⱼₖ, CartesianBoundary{d,LowerBoundary}}(),
+                    )
+                )
+            end
+        end
+    end
+
+    return c
+end
+
+"""
+    boundary_identifiers(a::CartesianAtlas)
+
+All non-connected boundaries of the charts of `a`.
+"""
+function boundary_identifiers(a::CartesianAtlas)
+    bs = MultiBlockBoundary[]
+
+    for d ∈ 1:ndims(charts(a))
+        Is = eachslice(CartesianIndices(charts(a)); dims=d)
+
+        for (i,b) ∈ ((1,LowerBoundary),(length(Is),UpperBoundary)) # For first and last slice
+            for jk ∈ eachindex(Is[i]) # For each block in slice
+                Iᵢⱼₖ = Tuple(Is[i][jk])
+                push!(bs,
+                    MultiBlockBoundary{Iᵢⱼₖ,   CartesianBoundary{d,b}}(),
+                )
+            end
+        end
+    end
+
+    return bs
+end
+
+
+"""
+    UnstructuredAtlas{C<:Chart, CN<:Tuple{MultiBlockBoundary,MultiBlockBoundary}, ...} <: Atlas
+
+An atlas with connections determined by a vector `MultiBlockBoundary` pairs.
+"""
+struct UnstructuredAtlas{C<:Chart, CN<:Tuple{MultiBlockBoundary,MultiBlockBoundary}, CV<:AbstractVector{C}, CNV<:AbstractVector{CN}} <: Atlas
+    charts::CV
+    connections::CNV
+end
+
+charts(a::UnstructuredAtlas) = a.charts
+connections(a::UnstructuredAtlas) = a.connections
+
+"""
+    boundary_identifiers(a::UnstructuredAtlas)
+
+All non-connected boundaries of the charts of `a`.
+"""
+function boundary_identifiers(a::UnstructuredAtlas)
+    bs = MultiBlockBoundary[]
+
+    for (i,c) ∈ enumerate(charts(a))
+        for b ∈ boundary_identifiers(c)
+            mbb = MultiBlockBoundary{i,typeof(b)}()
+
+            if !any(cn->mbb∈cn, connections(a))
+                push!(bs, mbb)
+            end
+        end
+    end
+
+    return bs
+end

--- a/test/Grids/equidistant_grid_test.jl	Tue Apr 29 08:48:58 2025 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -163,6 +163,15 @@
         @test equidistant_grid(unitinterval(),3) == equidistant_grid(0,1,3)
         @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η]
+        end
+        Grids.jacobian(c::typeof(c), ξ̄) = @SMatrix[2 0; 0 3]
+
+        @test equidistant_grid(c, 5, 4) isa Grid
+    end
 end
 
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Grids/manifolds_test.jl	Tue Apr 29 08:59:29 2025 +0200
@@ -0,0 +1,217 @@
+using Test
+
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
+using Diffinitive.LazyTensors
+
+using StaticArrays
+
+west = CartesianBoundary{1,LowerBoundary}
+east = CartesianBoundary{1,UpperBoundary}
+south = CartesianBoundary{2,LowerBoundary}
+north = CartesianBoundary{2,UpperBoundary}
+bottom = CartesianBoundary{3, LowerBoundary}
+top = CartesianBoundary{3, UpperBoundary}
+
+@testset "Chart" begin
+    X(ξ) = 2ξ
+    Grids.jacobian(::typeof(X), ξ) = @SMatrix[2 0; 0 2]
+    c = Chart(X, unitsquare())
+    @test c isa Chart{2}
+    @test parameterspace(c) == unitsquare()
+    @test ndims(c) == 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
+
+@testset "CartesianAtlas" begin
+    @testset "Constructors" begin
+        c = Chart(identity, unitsquare())
+        @test CartesianAtlas([c c; c c]) isa Atlas
+
+        c2 = Chart(x->2x, unitsquare())
+        @test CartesianAtlas([c c2; c2 c]) isa CartesianAtlas
+        @test CartesianAtlas(@SMatrix[c c; c c]) isa CartesianAtlas
+        @test CartesianAtlas(@SMatrix[c c2; c2 c]) isa CartesianAtlas
+    end
+
+    @testset "Getters" begin
+        c = Chart(identity, unitsquare())
+        a = CartesianAtlas([c c; c c])
+        @test charts(a) == [c c; c c]
+    end
+
+    @testset "connections" begin
+        # 2D
+        a = CartesianAtlas(fill(Chart(identity, unitsquare()), 2,3))
+
+        @test Set(connections(a)) == Set([
+            (MultiBlockBoundary{(1,1), east}(), MultiBlockBoundary{(2,1), west}()),
+            (MultiBlockBoundary{(1,1), north}(), MultiBlockBoundary{(1,2), south}()),
+            (MultiBlockBoundary{(2,1), north}(), MultiBlockBoundary{(2,2), south}()),
+            (MultiBlockBoundary{(1,2), east}(), MultiBlockBoundary{(2,2), west}()),
+            (MultiBlockBoundary{(1,2), north}(), MultiBlockBoundary{(1,3), south}()),
+            (MultiBlockBoundary{(2,2), north}(), MultiBlockBoundary{(2,3), south}()),
+            (MultiBlockBoundary{(1,3), east}(), MultiBlockBoundary{(2,3), west}()),
+        ])
+
+        # 3D
+        a = CartesianAtlas(fill(Chart(identity, unitcube()), 2,2,3))
+        @test Set(connections(a)) == Set([
+            (MultiBlockBoundary{(1,1,1), east}(),  MultiBlockBoundary{(2,1,1), west}()),
+            (MultiBlockBoundary{(1,1,1), north}(), MultiBlockBoundary{(1,2,1), south}()),
+            (MultiBlockBoundary{(2,1,1), north}(), MultiBlockBoundary{(2,2,1), south}()),
+            (MultiBlockBoundary{(1,2,1), east}(),  MultiBlockBoundary{(2,2,1), west}()),
+
+            (MultiBlockBoundary{(1,1,2), east}(),  MultiBlockBoundary{(2,1,2), west}()),
+            (MultiBlockBoundary{(1,1,2), north}(), MultiBlockBoundary{(1,2,2), south}()),
+            (MultiBlockBoundary{(2,1,2), north}(), MultiBlockBoundary{(2,2,2), south}()),
+            (MultiBlockBoundary{(1,2,2), east}(),  MultiBlockBoundary{(2,2,2), west}()),
+
+            (MultiBlockBoundary{(1,1,3), east}(),  MultiBlockBoundary{(2,1,3), west}()),
+            (MultiBlockBoundary{(1,1,3), north}(), MultiBlockBoundary{(1,2,3), south}()),
+            (MultiBlockBoundary{(2,1,3), north}(), MultiBlockBoundary{(2,2,3), south}()),
+            (MultiBlockBoundary{(1,2,3), east}(),  MultiBlockBoundary{(2,2,3), west}()),
+
+            (MultiBlockBoundary{(1,1,1), top}(), MultiBlockBoundary{(1,1,2), bottom}()),
+            (MultiBlockBoundary{(2,1,1), top}(), MultiBlockBoundary{(2,1,2), bottom}()),
+            (MultiBlockBoundary{(1,2,1), top}(), MultiBlockBoundary{(1,2,2), bottom}()),
+            (MultiBlockBoundary{(2,2,1), top}(), MultiBlockBoundary{(2,2,2), bottom}()),
+
+            (MultiBlockBoundary{(1,1,2), top}(), MultiBlockBoundary{(1,1,3), bottom}()),
+            (MultiBlockBoundary{(2,1,2), top}(), MultiBlockBoundary{(2,1,3), bottom}()),
+            (MultiBlockBoundary{(1,2,2), top}(), MultiBlockBoundary{(1,2,3), bottom}()),
+            (MultiBlockBoundary{(2,2,2), top}(), MultiBlockBoundary{(2,2,3), bottom}()),
+        ])
+    end
+
+    @testset "boundary_identifiers" begin
+        # 2D
+        a = CartesianAtlas(fill(Chart(identity, unitcube()), 2,3))
+        @test Set(boundary_identifiers(a)) == Set([
+            MultiBlockBoundary{(1,1), south}(),
+            MultiBlockBoundary{(2,1), south}(),
+            MultiBlockBoundary{(2,1), east}(),
+            MultiBlockBoundary{(2,2), east}(),
+            MultiBlockBoundary{(2,3), east}(),
+            MultiBlockBoundary{(1,3), north}(),
+            MultiBlockBoundary{(2,3), north}(),
+            MultiBlockBoundary{(1,1), west}(),
+            MultiBlockBoundary{(1,2), west}(),
+            MultiBlockBoundary{(1,3), west}(),
+        ])
+
+        # 3D
+        a = CartesianAtlas(fill(Chart(identity, unitsquare()), 2,2,3))
+        @test Set(boundary_identifiers(a)) == Set([
+            MultiBlockBoundary{(1,1,1), bottom}(),
+            MultiBlockBoundary{(2,1,1), bottom}(),
+            MultiBlockBoundary{(1,2,1), bottom}(),
+            MultiBlockBoundary{(2,2,1), bottom}(),
+
+            MultiBlockBoundary{(1,1,3), top}(),
+            MultiBlockBoundary{(2,1,3), top}(),
+            MultiBlockBoundary{(1,2,3), top}(),
+            MultiBlockBoundary{(2,2,3), top}(),
+
+            MultiBlockBoundary{(1,1,1), west}(),
+            MultiBlockBoundary{(1,2,1), west}(),
+            MultiBlockBoundary{(1,1,2), west}(),
+            MultiBlockBoundary{(1,2,2), west}(),
+            MultiBlockBoundary{(1,1,3), west}(),
+            MultiBlockBoundary{(1,2,3), west}(),
+
+            MultiBlockBoundary{(2,1,1), east}(),
+            MultiBlockBoundary{(2,2,1), east}(),
+            MultiBlockBoundary{(2,1,2), east}(),
+            MultiBlockBoundary{(2,2,2), east}(),
+            MultiBlockBoundary{(2,1,3), east}(),
+            MultiBlockBoundary{(2,2,3), east}(),
+
+            MultiBlockBoundary{(1,1,1), south}(),
+            MultiBlockBoundary{(2,1,1), south}(),
+            MultiBlockBoundary{(1,1,2), south}(),
+            MultiBlockBoundary{(2,1,2), south}(),
+            MultiBlockBoundary{(1,1,3), south}(),
+            MultiBlockBoundary{(2,1,3), south}(),
+
+            MultiBlockBoundary{(1,2,1), north}(),
+            MultiBlockBoundary{(2,2,1), north}(),
+            MultiBlockBoundary{(1,2,2), north}(),
+            MultiBlockBoundary{(2,2,2), north}(),
+            MultiBlockBoundary{(1,2,3), north}(),
+            MultiBlockBoundary{(2,2,3), north}(),
+        ])
+    end
+end
+
+@testset "UnstructuredAtlas" begin
+    @testset "Constructors" begin
+        c1 = Chart(identity, unitsquare())
+        c2 = Chart(x->2x, unitsquare())
+        cn = [
+            (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()),
+            (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()),
+            (MultiBlockBoundary{2, north}(),  MultiBlockBoundary{3, south}()),
+        ]
+
+        @test UnstructuredAtlas([c1, c1, c1], cn) isa UnstructuredAtlas
+        @test UnstructuredAtlas([c1, c2, c1, c2], cn) isa UnstructuredAtlas
+
+
+        cn = @SVector[
+            (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()),
+            (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()),
+            (MultiBlockBoundary{2, north}(),  MultiBlockBoundary{3, south}()),
+        ]
+        @test UnstructuredAtlas(@SVector[c1, c1, c1], cn) isa UnstructuredAtlas
+        @test UnstructuredAtlas(@SVector[c1, c2, c1, c2], cn) isa UnstructuredAtlas
+    end
+
+    @testset "Getters" begin
+        c = Chart(identity, unitsquare())
+        cn = [
+            (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()),
+            (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()),
+            (MultiBlockBoundary{2, north}(),  MultiBlockBoundary{3, south}()),
+        ]
+
+        a = UnstructuredAtlas([c, c, c], cn)
+
+        @test charts(a) == [c,c,c]
+        @test connections(a) == cn
+    end
+
+    @testset "boundary_identifiers" begin
+        c = Chart(identity, unitsquare())
+        cn = [
+            (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()),
+            (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()),
+            (MultiBlockBoundary{2, north}(),  MultiBlockBoundary{3, south}()),
+        ]
+
+        a = UnstructuredAtlas([c, c, c], cn)
+
+        @test Set(boundary_identifiers(a)) == Set([
+            MultiBlockBoundary{1, west}(),
+            MultiBlockBoundary{1, south}(),
+            MultiBlockBoundary{2, south}(),
+            MultiBlockBoundary{2, east}(),
+            MultiBlockBoundary{3, north}(),
+            MultiBlockBoundary{3, east}(),
+        ])
+
+    end
+end