--- a/src/Grids/Grids.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/src/Grids/Grids.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -20,6 +20,15 @@
 export unitcube
 export unithyperbox
 
+export verticies
+export unittriangle
+export unittetrahedron
+export unitsimplex
+
+export Chart
+export ConcreteChart
+export parameterspace
+
 # Grid
 export Grid
 export coordinate_size

--- a/src/Grids/equidistant_grid.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -127,6 +127,10 @@
 	return EquidistantGrid(range(limit_lower, limit_upper, length=size)) # TBD: Should it use LinRange instead?
 end
 
+
+equidistant_grid(hb::HyperBox, dims::Vararg{Int}) = equidistant_grid(hb.a, hb.b, dims...)
+
+
 CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
 
 

--- a/src/Grids/manifolds.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/src/Grids/manifolds.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -18,6 +18,7 @@
 [`Simplex`](@ref),
 """
 abstract type ParameterSpace{D} end
+Base.ndims(::ParameterSpace{D}) where D = D
 
 struct HyperBox{T,D} <: ParameterSpace{D}
     a::SVector{D,T}
@@ -43,53 +44,60 @@
 unithyperbox(D) = unithyperbox(Float64,D)
 
 
-struct Simplex{T,D} <: ParameterSpace{D}
-    verticies::NTuple{D,SVector{D,T}}
+struct Simplex{T,D,NV} <: ParameterSpace{D}
+    verticies::NTuple{NV,SVector{D,T}}
 end
 
+Simplex(verticies::Vararg{AbstractArray}) = Simplex(Tuple(SVector(v...) for v ∈ verticies))
+
+verticies(s::Simplex) = s.verticies
+
 Triangle{T} = Simplex{T,2}
 Tetrahedron{T} = Simplex{T,3}
 
-unittriangle(T) = unitsimplex(T,2)
-unittetrahedron(T) = unitsimplex(T,3)
+unittriangle(T=Float64) = unitsimplex(T,2)
+unittetrahedron(T=Float64) = unitsimplex(T,3)
 function unitsimplex(T,D)
     z = @SVector zeros(T,D)
     unitelement = one(eltype(z))
-    verticies = ntuple(i->setindex(z, unitelement, i), 4)
-    return Simplex(verticies)
+    verticies = ntuple(i->setindex(z, unitelement, i), D)
+    return Simplex((z,verticies...))
 end
-
+unitsimplex(D) = unitsimplex(Float64, D)
 
 """
+    Chart{D}
 
 A parametrized description of a manifold or part of a manifold.
-
-Should implement a methods for
-* `parameterspace`
-* `(::Chart)(ξs...)`
-
-There is a default implementation for `(::Chart{D})(::SVector{D})`
 """
-abstract type Chart{D} end
-# abstract type Chart{D,R} end
+struct Chart{D, PST<:ParameterSpace{D}, MT}
+    mapping::MT
+    parameterspace::PST
+end
 
 domain_dim(::Chart{D}) where D = D
-# range_dim(::Chart{D,R}) where {D,R} = R
+(c::Chart)(ξ) = c.mapping(ξ)
+parameterspace(c::Chart) = c.parameterspace
 
 """
-The parameterspace of a chart
-"""
-function parameterspace end
-
-(c::Chart{D})(x̄::SVector{D}) where D = c(x̄...)
-
+    jacobian(c::Chart, ξ)
 
-struct ConcereteChart{PST<:ParameterSpace, MT}
-    parameterspace::PST
-    mapping::MT
-end
-
-(c::Chart)(x̄) = c.mapping(x̄)
+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,ξ)`.
+"""
+jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ)
 
 
 """

--- a/src/Grids/mapped_grid.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -53,7 +53,7 @@
 
 function mapped_grid(x, J, size...)
     D = length(size)
-    lg = equidistant_grid(size, ntuple(i->0., D), ntuple(i->1., D))
+    lg = equidistant_grid(ntuple(i->0., D), ntuple(i->1., D), size...)
     return MappedGrid(
         lg,
         map(x,lg),
@@ -61,6 +61,15 @@
     )
 end
 
+function mapped_grid(c::Chart, size...)
+    lg = equidistant_grid(parameterspace(c), size...)
+    return MappedGrid(
+        lg,
+        map(c,lg),
+        map(ξ->jacobian(c, ξ), lg),
+    )
+end
+
 function jacobian_determinant(g::MappedGrid)
     return map(jacobian(g)) do ∂x∂ξ
         det(∂x∂ξ)

--- a/test/Grids/equidistant_grid_test.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -2,6 +2,7 @@
 using Test
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
+using StaticArrays
 
 
 @testset "EquidistantGrid" begin
@@ -143,6 +144,13 @@
             @test [gp[i]...] ≈ [p[i]...] atol=5e-13
         end
     end
+
+
+    @testset "equidistant_grid(::ParameterSpace)" begin
+        ps = HyperBox((0,0),(2,1))
+
+        @test equidistant_grid(ps, 3,4) == equidistant_grid((0,0), (2,1), 3,4)
+    end
 end
 
 

--- a/test/Grids/manifolds_test.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/test/Grids/manifolds_test.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -6,6 +6,11 @@
 
 # using StaticArrays
 
+@testset "ParameterSpace" begin
+    @test ndims(HyperBox([1,1], [2,2])) == 2
+    @test ndims(unittetrahedron()) == 3
+end
+
 @testset "HyperBox" begin
     @test HyperBox([1,1], [2,2]) isa HyperBox{Int, 2}
 
@@ -31,11 +36,24 @@
 
 @testset "Simplex" begin
     @test Simplex([1,2], [3,4]) isa Simplex{Int, 2}
-    @test Simplex([1,2,3], [4,5,6]) isa Simplex{Int, 3}
+    @test Simplex([1,2,3], [4,5,6],[1,1,1]) isa Simplex{Int, 3}
+
+    @test verticies(Simplex([1,2], [3,4])) == ([1,2], [3,4])
+
+    @test unittriangle() isa Simplex{Float64,2}
+    @test verticies(unittriangle()) == ([0,0], [1,0], [0,1])
+
+    @test unittetrahedron() isa  Simplex{Float64,3}
+    @test verticies(unittetrahedron()) == ([0,0,0], [1,0,0], [0,1,0],[0,0,1])
+
+    @test unitsimplex(4) isa Simplex{Float64,4}
 end
 
 @testset "Chart" begin
-
+    c = Chart(x->2x, unitsquare())
+    @test c isa Chart{2}
+    @test c([3,2]) == [6,4]
+    @test parameterspace(c) == unitsquare()
 end
 
 @testset "Atlas" begin

--- a/test/Grids/mapped_grid_test.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/test/Grids/mapped_grid_test.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -4,7 +4,7 @@
 using StaticArrays
 
 @testset "MappedGrid" begin
-    lg = equidistant_grid((11,11), (0,0), (1,1)) # TODO: Change dims of the grid to be different
+    lg = equidistant_grid((0,0), (1,1), 11, 11) # TODO: Change dims of the grid to be different
     x̄ = map(ξ̄ -> 2ξ̄, lg)
     J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
     mg = MappedGrid(lg, x̄, J)
@@ -64,7 +64,7 @@
 
     @testset "Iterator interface" begin
         sg = MappedGrid(
-            equidistant_grid((15,11), (0,0), (1,1)),
+            equidistant_grid((0,0), (1,1), 15, 11),
             map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11)
         )
 
@@ -180,7 +180,18 @@
     mg = mapped_grid(x̄, J, 10, 11)
     @test mg isa MappedGrid{SVector{2,Float64}, 2}
 
-    lg = equidistant_grid((10,11), (0,0), (1,1))
+    lg = equidistant_grid((0,0), (1,1), 10, 11)
     @test logicalgrid(mg) == lg
     @test collect(mg) == map(x̄, lg)
+
+
+    @testset "mapped_grid(::Chart)" begin
+        c = Chart(unitsquare()) do (ξ,η)
+            @SVector[2ξ, 3η]
+        end
+        Grids.jacobian(c::typeof(c), ξ̄) = @SMatrix[2 0; 0 3]
+
+        @test mapped_grid(c, 5, 4) isa Grid
+    end
+
 end