Mercurial > repos > public > sbplib_julia

changeset 1785:96f8cad255b4 feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 16 Sep 2024 10:33:47 +0200
parents f3d7e2d7a43f (current diff) c5070edd0ebb (diff)
children 1f42944d4a72
files Manifest.toml Project.toml src/Grids/Grids.jl
diffstat 22 files changed, 136 insertions(+), 824 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Wed Sep 11 16:26:19 2024 +0200
+++ b/Manifest.toml	Mon Sep 16 10:33:47 2024 +0200
@@ -2,66 +2,11 @@
 
 julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "32fac879810099260f177c27318d3f26de4a00cc"
-
-[[deps.Adapt]]
-deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
-uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.4"
-weakdeps = ["StaticArrays"]
-
-    [deps.Adapt.extensions]
-    AdaptStaticArraysExt = "StaticArrays"
-
-[[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76"
-uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.16.0"
-
-    [deps.ArrayInterface.extensions]
-    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
-    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
-    ArrayInterfaceCUDAExt = "CUDA"
-    ArrayInterfaceCUDSSExt = "CUDSS"
-    ArrayInterfaceChainRulesExt = "ChainRules"
-    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
-    ArrayInterfaceReverseDiffExt = "ReverseDiff"
-    ArrayInterfaceSparseArraysExt = "SparseArrays"
-    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
-    ArrayInterfaceTrackerExt = "Tracker"
-
-    [deps.ArrayInterface.weakdeps]
-    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
-    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
-    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
-    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
-    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+project_hash = "07350208c6e9bd0ec3979df9ac99bb401ac56208"
 
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
-[[deps.CommonWorldInvalidations]]
-git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0"
-uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8"
-version = "1.0.0"
-
-[[deps.Compat]]
-deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.16.0"
-weakdeps = ["Dates", "LinearAlgebra"]
-
-    [deps.Compat.extensions]
-    CompatLinearAlgebraExt = "LinearAlgebra"
-
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
@@ -71,11 +16,6 @@
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
-[[deps.IfElse]]
-git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
-uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
-version = "0.1.1"
-
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
@@ -83,15 +23,6 @@
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[deps.OffsetArrays]]
-git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
-uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.14.1"
-weakdeps = ["Adapt"]
-
-    [deps.OffsetArrays.extensions]
-    OffsetArraysAdaptExt = "Adapt"
-
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
@@ -117,33 +48,10 @@
 deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[[deps.Requires]]
-deps = ["UUIDs"]
-git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
-uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.3.0"
-
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
 
-[[deps.Static]]
-deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"]
-git-tree-sha1 = "87d51a3ee9a4b0d2fe054bdd3fc2436258db2603"
-uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "1.1.1"
-
-[[deps.StaticArrayInterface]]
-deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Static"]
-git-tree-sha1 = "96381d50f1ce85f2663584c8e886a6ca97e60554"
-uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718"
-version = "1.8.0"
-weakdeps = ["OffsetArrays", "StaticArrays"]
-
-    [deps.StaticArrayInterface.extensions]
-    StaticArrayInterfaceOffsetArraysExt = "OffsetArrays"
-    StaticArrayInterfaceStaticArraysExt = "StaticArrays"
-
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
 git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50"
@@ -168,16 +76,6 @@
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 version = "1.0.3"
 
-[[deps.TiledIteration]]
-deps = ["OffsetArrays", "StaticArrayInterface"]
-git-tree-sha1 = "1176cc31e867217b06928e2f140c90bd1bc88283"
-uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
-version = "0.5.0"
-
-[[deps.UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
 [[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
--- a/Notes.md	Wed Sep 11 16:26:19 2024 +0200
+++ b/Notes.md	Mon Sep 16 10:33:47 2024 +0200
@@ -118,19 +118,10 @@
  Maybe if we should have dynamic sizing it could be only for the range. `domain_size` would not be implemented. And the `range_size` would be a function of a vector that the LazyTensor is applied to.
 
 ## Reasearch and thinking
- - [ ] Use a trait to indicate that a LazyTensor har the same range and domain?
  - [ ] Check how the native julia doc generator works
     - [ ] Check if Vidars design docs fit in there
  - [ ] Create a macro @lazy which replaces a binary op (+,-) by its lazy equivalent? Would be a neat way to indicate which evaluations are lazy without cluttering/confusing with special characters.
- - [ ] Dispatch on Lower() instead of the type Lower so `::Lower` instead of `::Type{Lower}` ???
- 	Seems better unless there is some specific reason to use the type instead of the value.
  - [ ] Can we have a trait to tell if a LazyTensor is transposable?
- - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input?
- - [ ] Figure out how to treat the borrowing parameters of operators. Include in into the struct? Expose via function dispatched on the operator type and grid?
-
-## Identifiers for regions
-The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probably be included in the grid module. This allows new grid types to come with their own regions.
-We implement this by refactoring RegionIndices to be agnostic to the region types and then moving the actual types to Grids.
 
 ## Regions and tensormappings
 - [ ] Use a trait to indicate if a LazyTensor uses indices with regions.
@@ -182,15 +173,14 @@
 Preferably dimensions and sizes should be checked when lazy objects are created, for example TensorApplication, TensorComposition and so on. If dimension checks decreases performance we can make them skippable later.
 
 ## Changes to `eval_on`
-There are reasons to replace `eval_on` with regular `map` from Base, and implement a kind of lazy map perhaps `lmap` that work on indexable collections.
-
-The benefit of doing this is that we can treat grids as gridfunctions for the coordinate function, and get a more flexible tool. For example `map`/`lmap` can then be used both to evaluate a function on the grid but also get a component of a vector valued grid function or similar.
+There are reasons to replace `eval_on` with regular `map` from Base, and
+implement a kind of lazy map perhaps `lmap` that work on indexable
+collections.
 
-A question is how and if we should implement `map`/`lmap` for functions like `(x,y)->x*y` or stick to just using vector inputs. There are a few options.
-
-* use `Base.splat((x,y)->x*y)` with the single argument `map`/`lmap`.
-* implement a kind of `unzip` function to get iterators for each component, which can then be used with the multiple-iterators-version of `map`/`lmap`.
-* Inspect the function in the `map`/`lmap` function to determine which matches.
+The benefit of doing this is that we can treat grids as gridfunctions for the
+coordinate function, and get a more flexible tool. For example `map`/`lmap`
+can then be used both to evaluate a function on the grid but also get a
+component of a vector valued grid function or similar.
 
 Below is a partial implementation of `lmap` with some ideas
 ```julia
@@ -220,7 +210,10 @@
 lmap(f,  I) = LazyIndexableMap(f,I)
 ```
 
-The interaction of the map methods with the probable design of multiblock functions involving nested indecies complicate the picture slightly. It's clear at the time of writing how this would work with `Base.map`. Perhaps we want to implement our own versions of both eager and lazy map.
+The interaction of the map methods with the probable design of multiblock
+functions involving nested indecies complicate the picture slightly. It's
+unclear at the time of writing how this would work with `Base.map`. Perhaps we
+want to implement our own versions of both eager and lazy map.
 
 
 ### 2024-04
@@ -291,16 +284,6 @@
 Grid-funktioner har typen `AbstractArray{T,2} where T`.
 `T` kan vara lite vad som helst, tillexemel en SVector eller Array, eller Tuple. Tensoroperatorerna bryr sig inte om exakt vad det är, mer än att typen måste stödja de operationer som operatorn använder.
 
-En nackdel kan vara hur man ska få ut gridfunktionen för tex andra komponenten.
-
-Syntax:
-```
-f(x̄) = x̄
-gf = evalOn(g, f)
-gf[2,3] # x̄ för en viss gridpunkt
-gf[2,3][2] # x̄[2] för en viss gridpunkt
-```
-
 ### Tensor operatorer
 Vi kan ha tensor-operatorer som agerar på ett skalärt fält och ger ett vektorfält eller tensorfält.
 Vi kan också ha tensor-operatorer som agerar på ett vektorfält eller tensorfält och ger ett skalärt fält.
@@ -316,16 +299,6 @@
 
 Kanske kan man implementera `⋅(tm::LazyTensor{R,D}, v::AbstractArray{T,D})` där T är en AbstractArray, tm på något sätt har komponenter, lika många som T har element.
 
-### Komponenter som gridfunktioner
-En viktig operation för vektorfält är att kunna få ut komponenter som grid-funktioner. Detta behöver antagligen kunna ske lazy.
-Det finns ett par olika lösningar:
-* Använda map eller en lazy map (se diskussion om eval_on)
-* Implementera en egen typ av view som tar hand om detta. Eller Accessors.jl?
-* Använda en LazyTensor
-* Någon typ av lazy-broadcast
-* En lazy array som applicerar en funktion för varje element.
-
-
 ### Prestanda-aspekter
 [Vidar, Discord, 2023-03-03]
 Typiskt sett finns det två sätt att representera vektorvärda gridfunktioner AbstractArray{T,Dim} där T är en vektor över komponenterna. Man skulle alltså i 1D ha
@@ -392,3 +365,38 @@
 stencil application on the bugfix/sbp_operators/stencil_return_type branch
 there seemed to be some strange results where such errors could be the
 culprit.
+
+
+## Tiled loops and regions in apply
+There should be easy ways to use functionalty splitting the application of a lazy array into regions and using tiled iteration. This could make the application more efficient by reducing branching and improving cache usage in the tight loop. On commit f215ac2a5c66 and before there were some early tests regarding this in a DiffOp submodule.
+
+The main ideas were:
+```julia
+function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
+    apply_region!(D, u, v, Lower, Lower)
+    apply_region!(D, u, v, Lower, Interior)
+    apply_region!(D, u, v, Lower, Upper)
+    apply_region!(D, u, v, Interior, Lower)
+    apply_region!(D, u, v, Interior, Interior)
+    apply_region!(D, u, v, Interior, Upper)
+    apply_region!(D, u, v, Upper, Lower)
+    apply_region!(D, u, v, Upper, Interior)
+    apply_region!(D, u, v, Upper, Upper)
+    return nothing
+end
+```
+
+```julia
+using TiledIteration
+function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
+    ri = regionindices(D.grid.size, closuresize(D.op), (r1,r2))
+    # TODO: Pass Tilesize to function
+    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
+        for j ∈ tileaxs[2], i ∈ tileaxs[1]
+            I = ri[i,j]
+            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
+        end
+    end
+    return nothing
+end
+```
--- a/Project.toml	Wed Sep 11 16:26:19 2024 +0200
+++ b/Project.toml	Mon Sep 16 10:33:47 2024 +0200
@@ -7,7 +7,6 @@
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
-TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
--- a/TODO.md	Wed Sep 11 16:26:19 2024 +0200
+++ b/TODO.md	Mon Sep 16 10:33:47 2024 +0200
@@ -4,11 +4,7 @@
  - [ ] Split up Notes.md in several files
 
 ## Coding
- - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
- - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?)
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
- - [ ] Use `@inferred` in a lot of tests.
- - [ ] Replace `@inferred` tests with a benchmark suite that automatically tests for regressions.
  - [ ] Make sure we are setting tolerances in tests in a consistent way
  - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index
  - [ ] Clean up RegionIndices
@@ -16,7 +12,6 @@
     2. [ ] Update RegionIndices accordingly
     3. [ ] Fix the rest of the library
     Should getregion also work for getregion(::Colon,...)
- - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`?
  - [ ] Add custom pretty printing to LazyTensors/SbpOperators to enhance readability of e.g error messages.
        See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing)
  - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe
--- a/TimeStepper.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,60 +0,0 @@
-abstract type TimeStepper end
-
-# Returns v and t
-function getState(ts::TimeStepper)
-	error("not implemented")
-end
-
-
-function step!(ts::TimeStepper)
-	error("not implemented")
-end
-
-function stepN(ts::TimeStepper,N::Int)
-	for i ∈ 1:N
-		ts.step()
-	end
-end
-
-function stepTo(ts::TimeStepper)
-	error("Not yet implemented")
-end
-
-function evolve(ts::TimeStepper)
-	error("Not yet implemented")
-end
-
-
-mutable struct Rk4 <: TimeStepper
-	F::Function
-	k::Real
-	v::Vector
-	t::Real
-	n::UInt
-
-	function Rk4(F::Function,k::Real,v0::Vector,t0::Real)
-		# TODO: Check that F has two inputs and one output
-		v = v0
-		t = t0
-		n = 0
-		return new(F,k,v,t,n)
-	end
-end
-
-function getState(ts::Rk4)
-	return ts.t, ts.v
-end
-
-function step!(ts::Rk4)
-    k1 = ts.F(ts.v,ts.t)
-	k2 = ts.F(ts.v+0.5*ts.k*k1,ts.t+0.5*ts.k)
-	k3 = ts.F(ts.v+0.5*ts.k*k2,ts.t+0.5*ts.k)
-    k4 = ts.F(ts.v+    ts.k*k3,ts.t+    ts.k)
-    ts.v  = ts.v + (1/6)*(k1+2*(k2+k3)+k4)*ts.k
-
-	ts.n = ts.n + 1
-	ts.t = ts.t + ts.k
-
-	return nothing
-end
-
--- a/docs/make.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/docs/make.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -1,12 +1,10 @@
 using Documenter
 using Diffinitive
 
-using Diffinitive.DiffOps
 using Diffinitive.Grids
 using Diffinitive.LazyTensors
 using Diffinitive.RegionIndices
 using Diffinitive.SbpOperators
-using Diffinitive.StaticDicts
 
 sitename = "Diffinitive.jl"
 
@@ -34,11 +32,9 @@
     "matrix_and_tensor_representations.md",
     "Submodules" => [
         "submodules/grids.md",
-        "submodules/diff_ops.md",
         "submodules/lazy_tensors.md",
         "submodules/region_indices.md",
         "submodules/sbp_operators.md",
-        "submodules/static_dicts.md",
     ],
     "doc_index.md",
 ]
--- a/docs/src/submodules/diff_ops.md	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,23 +0,0 @@
-# DiffOps
-
-## Contents
-```@contents
-Pages = ["diff_ops.md"]
-```
-
-## Index
-```@index
-Pages = ["diff_ops.md"]
-```
-
-## Public interface
-```@autodocs
-Modules = [DiffOps]
-Private = false # Hide unexported objects
-```
-
-## Internal interface
-```@autodocs
-Modules = [DiffOps]
-Public = false # Hide exported objects
-```
--- a/docs/src/submodules/static_dicts.md	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,23 +0,0 @@
-# StaticDicts
-
-## Contents
-```@contents
-Pages = ["static_dicts.md"]
-```
-
-## Index
-```@index
-Pages = ["static_dicts.md"]
-```
-
-## Public interface
-```@autodocs
-Modules = [StaticDicts]
-Private = false # Hide unexported objects
-```
-
-## Internal interface
-```@autodocs
-Modules = [StaticDicts]
-Public = false # Hide exported objects
-```
--- a/plotDerivative.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,12 +0,0 @@
-g = sbp.Grid.EquidistantGrid((200,), (0.0,), (2pi,))
-op =sbp.readOperator("d2_4th.txt","h_4th.txt")
-Laplace = sbp.Laplace(g,1.0,op)
-
-init(x) = cos(x)
-v = sbp.Grid.evalOn(g,init)
-u = zeros(length(v))
-
-sbp.apply!(Laplace,u,v)
-
-@show u
-sbp.Grid.plotgridfunction(g,u)
--- a/plotDerivative2d.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,16 +0,0 @@
-include("sbpPlot.jl")
-
-g = sbp.Grid.EquidistantGrid((100,75), (0.0, 0.0), (2pi, 3/2*pi))
-op = sbp.readOperator("d2_4th.txt","h_4th.txt")
-Laplace = sbp.Laplace(g, 1.0, op)
-
-init(x,y) = sin(x) + cos(y)
-v = sbp.Grid.evalOn(g,init)
-u = zero(v)
-
-sbp.apply!(Laplace,u,v)
-
-#@show u
-#@show u'*u
-
-plotgridfunction(g,u)
--- a/src/DiffOps/DiffOps.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,103 +0,0 @@
-module DiffOps
-
-using Diffinitive.RegionIndices
-using Diffinitive.SbpOperators
-using Diffinitive.Grids
-using Diffinitive.LazyTensors
-
-"""
-    DiffOp
-
-Supertype of differential operator discretisations.
-The action of the DiffOp is defined in the method
-    apply(D::DiffOp, v::AbstractVector, I...)
-"""
-abstract type DiffOp end
-
-function apply end
-
-function matrixRepresentation(D::DiffOp)
-    error("not implemented")
-end
-
-abstract type DiffOpCartesian{Dim} <: DiffOp end
-
-# DiffOp must have a grid of dimension Dim!!!
-function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
-    for I ∈ eachindex(D.grid)
-        u[I] = apply(D, v, I)
-    end
-
-    return nothing
-end
-export apply!
-
-function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
-    apply_region!(D, u, v, Lower, Lower)
-    apply_region!(D, u, v, Lower, Interior)
-    apply_region!(D, u, v, Lower, Upper)
-    apply_region!(D, u, v, Interior, Lower)
-    apply_region!(D, u, v, Interior, Interior)
-    apply_region!(D, u, v, Interior, Upper)
-    apply_region!(D, u, v, Upper, Lower)
-    apply_region!(D, u, v, Upper, Interior)
-    apply_region!(D, u, v, Upper, Upper)
-    return nothing
-end
-
-# Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
-function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    for I ∈ regionindices(D.grid.size, closuresize(D.op), (r1,r2))
-        @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
-        @inbounds u[I] = apply(D, v, indextuple)
-    end
-    return nothing
-end
-export apply_region!
-
-function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
-    apply_region_tiled!(D, u, v, Lower, Lower)
-    apply_region_tiled!(D, u, v, Lower, Interior)
-    apply_region_tiled!(D, u, v, Lower, Upper)
-    apply_region_tiled!(D, u, v, Interior, Lower)
-    apply_region_tiled!(D, u, v, Interior, Interior)
-    apply_region_tiled!(D, u, v, Interior, Upper)
-    apply_region_tiled!(D, u, v, Upper, Lower)
-    apply_region_tiled!(D, u, v, Upper, Interior)
-    apply_region_tiled!(D, u, v, Upper, Upper)
-    return nothing
-end
-
-using TiledIteration
-function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    ri = regionindices(D.grid.size, closuresize(D.op), (r1,r2))
-    # TODO: Pass Tilesize to function
-    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
-        for j ∈ tileaxs[2], i ∈ tileaxs[1]
-            I = ri[i,j]
-            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
-        end
-    end
-    return nothing
-end
-export apply_region_tiled!
-
-function apply(D::DiffOp, v::AbstractVector)::AbstractVector
-    u = zeros(eltype(v), size(v))
-    apply!(D,v,u)
-    return u
-end
-
-# TODO: This conflicts with LazyTensors. Shouldn't DiffOps be LazyTensorOperators and use that apply?
-# export apply
-
-
-"""
-    BoundaryCondition
-A BoundaryCondition should implement the method
-    sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
-"""
-abstract type BoundaryCondition end
-
-
-end # module
--- a/src/Diffinitive.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/src/Diffinitive.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -1,16 +1,13 @@
 module Diffinitive
 
-include("StaticDicts/StaticDicts.jl")
 include("RegionIndices/RegionIndices.jl")
 include("LazyTensors/LazyTensors.jl")
 include("Grids/Grids.jl")
 include("SbpOperators/SbpOperators.jl")
-include("DiffOps/DiffOps.jl")
 
 export RegionIndices
 export LazyTensors
 export Grids
 export SbpOperators
-export DiffOps
 
 end
--- a/src/Grids/Grids.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/src/Grids/Grids.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -1,4 +1,5 @@
 # TODO: Double check that the interfaces for indexing and iterating are fully implemented and tested for all grids.
+# Review: Address this todo?
 module Grids
 
 using Diffinitive.LazyTensors
@@ -64,9 +65,8 @@
 # MappedGrid
 export MappedGrid
 export jacobian
-export logicalgrid
+export logical_grid
 export mapped_grid
-export jacobian_determinant
 export metric_tensor
 export metric_tensor_inverse
 
--- a/src/Grids/equidistant_grid.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -151,9 +151,8 @@
 	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...)
-# TODO: One dimensional grids shouldn't have vector eltype right?, Change here or in HyperBox?
+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})
     lg = equidistant_grid(parameterspace(c), dims...)
--- a/src/Grids/manifolds.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/src/Grids/manifolds.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -19,7 +19,21 @@
 """
 abstract type ParameterSpace{D} end
 Base.ndims(::ParameterSpace{D}) where D = D
-# TBD:  Should implement domain_dim?
+
+struct Interval{T} <: ParameterSpace{1}
+    a::T
+    b::T
+
+    function Interval(a,b)
+        a, b = promote(a, b)
+        new{typeof(a)}(a,b)
+    end
+end
+
+limits(i::Interval) = (i.a, i.b)
+
+unitinterval(T=Float64) = Interval(zero(T), one(T))
+
 
 struct HyperBox{T,D} <: ParameterSpace{D}
     a::SVector{D,T}
@@ -27,18 +41,17 @@
 end
 
 function HyperBox(a,b)
-    T = SVector{length(a)}
+    ET = promote_type(eltype(a),eltype(b))
+    T = SVector{length(a),ET}
     HyperBox(convert(T,a), convert(T,b))
 end
 
-Interval{T} = HyperBox{T,1}
 Rectangle{T} = HyperBox{T,2}
 Box{T} = HyperBox{T,3}
 
 limits(box::HyperBox, d) = (box.a[d], box.b[d])
 limits(box::HyperBox) = (box.a, box.b)
 
-unitinterval(T=Float64) = unithyperbox(T,1)
 unitsquare(T=Float64) = unithyperbox(T,2)
 unitcube(T=Float64) = unithyperbox(T,3)
 unithyperbox(T, D) = HyperBox((@SVector zeros(T,D)), (@SVector ones(T,D)))
@@ -49,7 +62,12 @@
     verticies::NTuple{NV,SVector{D,T}}
 end
 
-Simplex(verticies::Vararg{AbstractArray}) = Simplex(Tuple(SVector(v...) for v ∈ verticies))
+function Simplex(verticies::Vararg{AbstractArray})
+    ET = mapreduce(eltype,promote_type,verticies)
+    T = SVector{length(verticies[1]),ET}
+
+    return Simplex(Tuple(convert(T,v) for v ∈ verticies))
+end
 
 verticies(s::Simplex) = s.verticies
 
@@ -76,7 +94,7 @@
     parameterspace::PST
 end
 
-domain_dim(::Chart{D}) where D = D
+Base.ndims(::Chart{D}) where D = D
 (c::Chart)(ξ) = c.mapping(ξ)
 parameterspace(c::Chart) = c.parameterspace
 
@@ -131,6 +149,7 @@
 end
 
 charts(a::CartesianAtlas) = a.charts
+connections(a::CartesianAtlas) = nothing
 
 struct UnstructuredAtlas <: Atlas
     charts::Vector{Chart}
@@ -138,6 +157,7 @@
 end
 
 charts(a::UnstructuredAtlas) = a.charts
+connections(a::UnstructuredAtlas) = nothing
 
 
 ###
--- a/src/Grids/mapped_grid.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -5,20 +5,20 @@
 coordinates. The physical coordinates and the Jacobian are stored as grid
 functions corresponding to the logical grid.
 
-See also: [`logicalgrid`](@ref), [`jacobian`](@ref), [`metric_tensor`](@ref).
+See also: [`logical_grid`](@ref), [`jacobian`](@ref), [`metric_tensor`](@ref).
 """
-struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D}
-    logicalgrid::GT
+struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractMatrix{<:Any}, D}} <: Grid{T,D}
+    logical_grid::GT
     physicalcoordinates::CT
     jacobian::JT
 
     """
-        MappedGrid(logicalgrid, physicalcoordinates, jacobian)
+        MappedGrid(logical_grid, physicalcoordinates, jacobian)
 
     A MappedGrid with the given physical coordinates and jacobian.
     """
-    function MappedGrid(logicalgrid::GT, physicalcoordinates::CT, jacobian::JT) where {T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}}
-        if !(size(logicalgrid) == size(physicalcoordinates) == size(jacobian))
+    function MappedGrid(logical_grid::GT, physicalcoordinates::CT, jacobian::JT) where {T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractMatrix{<:Any}, D}}
+        if !(size(logical_grid) == size(physicalcoordinates) == size(jacobian))
             throw(ArgumentError("Sizes must match"))
         end
 
@@ -26,24 +26,24 @@
             throw(ArgumentError("The size of the jacobian must match the dimensions of the grid and coordinates"))
         end
 
-        return new{T,D,GT,CT,JT}(logicalgrid, physicalcoordinates, jacobian)
+        return new{T,D,GT,CT,JT}(logical_grid, physicalcoordinates, jacobian)
     end
 end
 
 function Base.:(==)(a::MappedGrid, b::MappedGrid)
-    same_logicalgrid = logicalgrid(a) == logicalgrid(b)
+    same_logical_grid = logical_grid(a) == logical_grid(b)
     same_coordinates = collect(a) == collect(b)
     same_jacobian = jacobian(a) == jacobian(b)
 
-    return same_logicalgrid && same_coordinates && same_jacobian
+    return same_logical_grid && same_coordinates && same_jacobian
 end
 
 """
-    logicalgrid(g::MappedGrid)
+    logical_grid(g::MappedGrid)
 
 The logical grid of `g`.
 """
-logicalgrid(g::MappedGrid) = g.logicalgrid
+logical_grid(g::MappedGrid) = g.logical_grid
 
 """
     jacobian(g::MappedGrid)
@@ -55,25 +55,27 @@
 
 # Indexing interface
 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...]
-Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid)
+Base.eachindex(g::MappedGrid) = eachindex(g.logical_grid)
 
-Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d)
-Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d)
+Base.firstindex(g::MappedGrid, d) = firstindex(g.logical_grid, d)
+Base.lastindex(g::MappedGrid, d) = lastindex(g.logical_grid, d)
 
 # Iteration interface
 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates)
 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state)
 
 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}()
-Base.length(g::MappedGrid) = length(g.logicalgrid)
-Base.size(g::MappedGrid) = size(g.logicalgrid)
-Base.size(g::MappedGrid, d) = size(g.logicalgrid, d)
+Base.length(g::MappedGrid) = length(g.logical_grid)
+Base.size(g::MappedGrid) = size(g.logical_grid)
+Base.size(g::MappedGrid, d) = size(g.logical_grid, d)
 
-boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logicalgrid)
-boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logicalgrid, id)
+boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logical_grid)
+boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logical_grid, id)
 
+# Review: Error when calling plot(boundary_grid(g, id))
+# Currently need to collect first, i.e., plot(collect(boundary_grid(g, id)))
 function boundary_grid(g::MappedGrid, id::TensorGridBoundary)
-    b_indices = boundary_indices(g.logicalgrid, id)
+    b_indices = boundary_indices(g.logical_grid, id)
 
     # Calculate indices of needed jacobian components
     D = ndims(g)
@@ -86,7 +88,7 @@
     boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...]
 
     return MappedGrid(
-        boundary_grid(g.logicalgrid, id),
+        boundary_grid(g.logical_grid, id),
         boundary_physicalcoordinates,
         boundary_jacobian,
     )
@@ -121,21 +123,6 @@
 end
 
 """
-    jacobian_determinant(g::MappedGrid)
-
-The jacobian determinant of `g` as a grid function.
-"""
-function jacobian_determinant(g::MappedGrid)
-    return map(jacobian(g)) do ∂x∂ξ
-        det(∂x∂ξ)
-    end
-end
-# TBD: Should this be changed to calculate sqrt(g) instead?
-#       This would make it well defined also for n-dim grids embedded in higher dimensions.
-# TBD: Is there a better name? metric_determinant?
-# TBD: Is the best option to delete it?
-
-"""
     metric_tensor(g::MappedGrid)
 
 The metric tensor of `g` as a grid function.
--- a/src/StaticDicts/StaticDicts.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,82 +0,0 @@
-module StaticDicts
-
-export StaticDict
-
-"""
-    StaticDict{K,V,N} <: AbstractDict{K,V}
-
-A static dictionary implementing the interface for an `AbstractDict`. A
-`StaticDict` is fully immutable and after creation no changes can be made.
-
-The immutable nature means that `StaticDict` can be compared with `===`, in
-constrast to regular `Dict` or `ImmutableDict` which can not. (See
-<https://github.com/JuliaLang/julia/issues/4648> for details.) One important
-aspect of this is that `StaticDict` can be used in a struct while still
-allowing the struct to be compared using the default implementation of `==` for
-structs.
-
-Lookups are done by linear search.
-
-Duplicate keys are not allowed and an error will be thrown if they are passed
-to the constructor.
-"""
-struct StaticDict{K,V,N} <: AbstractDict{K,V}
-    pairs::NTuple{N,Pair{K,V}}
-
-    function StaticDict{K,V}(pairs::Vararg{Pair,N}) where {K,V,N}
-        if !allunique(first.(pairs))
-            throw(DomainError(pairs, "keys must be unique"))
-        end
-        return new{K,V,N}(pairs)
-    end
-end
-
-function StaticDict(pairs::Vararg{Pair})
-    K = typejoin(firsttype.(pairs)...)
-    V = typejoin(secondtype.(pairs)...)
-    return StaticDict{K,V}(pairs...)
-end
-
-StaticDict(pairs::NTuple{N,Pair} where N) = StaticDict(pairs...)
-
-function Base.get(d::StaticDict, key, default)
-    for p ∈ d.pairs
-        if key == p.first
-            return p.second
-        end
-    end
-
-    return default
-end
-
-Base.iterate(d::StaticDict) = iterate(d.pairs)
-Base.iterate(d::StaticDict, state) = iterate(d.pairs,state)
-Base.length(d::StaticDict) = length(d.pairs)
-
-
-"""
-    merge(d1::StaticDict, d2::StaticDict)
-
-Merge two `StaticDict`. Repeating keys is considered and error. This may
-change in a future version.
-"""
-function Base.merge(d1::StaticDict, d2::StaticDict)
-    return StaticDict(d1.pairs..., d2.pairs...)
-end
-
-
-"""
-    firsttype(::Pair{T1,T2})
-
-The type of the first element in the pair.
-"""
-firsttype(::Pair{T1,T2}) where {T1,T2} = T1
-
-"""
-    secondtype(::Pair{T1,T2})
-
-The type of the secondtype element in the pair.
-"""
-secondtype(::Pair{T1,T2}) where {T1,T2}  = T2
-
-end # module
--- a/test/DiffOps/DiffOps_test.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,195 +0,0 @@
-using Test
-using Diffinitive.DiffOps
-using Diffinitive.Grids
-using Diffinitive.SbpOperators
-using Diffinitive.RegionIndices
-using Diffinitive.LazyTensors
-
-#
-# @testset "BoundaryValue" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
-#
-#     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
-#     e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
-#     e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
-#     e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = zeros(Float64, 4, 5)
-#     v[:,5] = [1, 2, 3,4]
-#     v[:,4] = [1, 2, 3,4]
-#     v[:,3] = [4, 5, 6, 7]
-#     v[:,2] = [7, 8, 9, 10]
-#     v[:,1] = [10, 11, 12, 13]
-#
-#     @test e_w  isa LazyTensor{T,2,1} where T
-#     @test e_w' isa LazyTensor{T,1,2} where T
-#
-#     @test domain_size(e_w, (3,2)) == (2,)
-#     @test domain_size(e_e, (3,2)) == (2,)
-#     @test domain_size(e_s, (3,2)) == (3,)
-#     @test domain_size(e_n, (3,2)) == (3,)
-#
-#     @test size(e_w'*v) == (5,)
-#     @test size(e_e'*v) == (5,)
-#     @test size(e_s'*v) == (4,)
-#     @test size(e_n'*v) == (4,)
-#
-#     @test collect(e_w'*v) == [10,7,4,1.0,1]
-#     @test collect(e_e'*v) == [13,10,7,4,4.0]
-#     @test collect(e_s'*v) == [10,11,12,13.0]
-#     @test collect(e_n'*v) == [1,2,3,4.0]
-#
-#     g_x = [1,2,3,4.0]
-#     g_y = [5,4,3,2,1.0]
-#
-#     G_w = zeros(Float64, (4,5))
-#     G_w[1,:] = g_y
-#
-#     G_e = zeros(Float64, (4,5))
-#     G_e[4,:] = g_y
-#
-#     G_s = zeros(Float64, (4,5))
-#     G_s[:,1] = g_x
-#
-#     G_n = zeros(Float64, (4,5))
-#     G_n[:,5] = g_x
-#
-#     @test size(e_w*g_y) == (UnknownDim,5)
-#     @test size(e_e*g_y) == (UnknownDim,5)
-#     @test size(e_s*g_x) == (4,UnknownDim)
-#     @test size(e_n*g_x) == (4,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken collect(e_w*g_y) == G_w
-#     @test_broken collect(e_e*g_y) == G_e
-#     @test_broken collect(e_s*g_x) == G_s
-#     @test_broken collect(e_n*g_x) == G_n
-# end
-#
-# @testset "NormalDerivative" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
-#
-#     d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
-#     d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
-#     d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
-#     d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
-#
-#
-#     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-#     v∂x = evalOn(g, (x,y)-> 2*x + y)
-#     v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
-#
-#     @test d_w  isa LazyTensor{T,2,1} where T
-#     @test d_w' isa LazyTensor{T,1,2} where T
-#
-#     @test domain_size(d_w, (3,2)) == (2,)
-#     @test domain_size(d_e, (3,2)) == (2,)
-#     @test domain_size(d_s, (3,2)) == (3,)
-#     @test domain_size(d_n, (3,2)) == (3,)
-#
-#     @test size(d_w'*v) == (6,)
-#     @test size(d_e'*v) == (6,)
-#     @test size(d_s'*v) == (5,)
-#     @test size(d_n'*v) == (5,)
-#
-#     @test collect(d_w'*v) ≈ v∂x[1,:]
-#     @test collect(d_e'*v) ≈ v∂x[5,:]
-#     @test collect(d_s'*v) ≈ v∂y[:,1]
-#     @test collect(d_n'*v) ≈ v∂y[:,6]
-#
-#
-#     d_x_l = zeros(Float64, 5)
-#     d_x_u = zeros(Float64, 5)
-#     for i ∈ eachindex(d_x_l)
-#         d_x_l[i] = op.dClosure[i-1]
-#         d_x_u[i] = -op.dClosure[length(d_x_u)-i]
-#     end
-#
-#     d_y_l = zeros(Float64, 6)
-#     d_y_u = zeros(Float64, 6)
-#     for i ∈ eachindex(d_y_l)
-#         d_y_l[i] = op.dClosure[i-1]
-#         d_y_u[i] = -op.dClosure[length(d_y_u)-i]
-#     end
-#
-#     function prod_matrix(x,y)
-#         G = zeros(Float64, length(x), length(y))
-#         for I ∈ CartesianIndices(G)
-#             G[I] = x[I[1]]*y[I[2]]
-#         end
-#
-#         return G
-#     end
-#
-#     g_x = [1,2,3,4.0,5]
-#     g_y = [5,4,3,2,1.0,11]
-#
-#     G_w = prod_matrix(d_x_l, g_y)
-#     G_e = prod_matrix(d_x_u, g_y)
-#     G_s = prod_matrix(g_x, d_y_l)
-#     G_n = prod_matrix(g_x, d_y_u)
-#
-#
-#     @test size(d_w*g_y) == (UnknownDim,6)
-#     @test size(d_e*g_y) == (UnknownDim,6)
-#     @test size(d_s*g_x) == (5,UnknownDim)
-#     @test size(d_n*g_x) == (5,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken collect(d_w*g_y) ≈ G_w
-#     @test_broken collect(d_e*g_y) ≈ G_e
-#     @test_broken collect(d_s*g_x) ≈ G_s
-#     @test_broken collect(d_n*g_x) ≈ G_n
-# end
-#
-# @testset "BoundaryQuadrature" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
-#
-#     H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
-#     H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
-#     H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
-#     H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-#
-#     function get_quadrature(N)
-#         qc = op.quadratureClosure
-#         q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
-#         @assert length(q) == N
-#         return q
-#     end
-#
-#     v_w = v[1,:]
-#     v_e = v[10,:]
-#     v_s = v[:,1]
-#     v_n = v[:,11]
-#
-#     q_x = spacing(g)[1].*get_quadrature(10)
-#     q_y = spacing(g)[2].*get_quadrature(11)
-#
-#     @test H_w isa TensorOperator{T,1} where T
-#
-#     @test domain_size(H_w, (3,)) == (3,)
-#     @test domain_size(H_n, (3,)) == (3,)
-#
-#     @test range_size(H_w, (3,)) == (3,)
-#     @test range_size(H_n, (3,)) == (3,)
-#
-#     @test size(H_w*v_w) == (11,)
-#     @test size(H_e*v_e) == (11,)
-#     @test size(H_s*v_s) == (10,)
-#     @test size(H_n*v_n) == (10,)
-#
-#     @test collect(H_w*v_w) ≈ q_y.*v_w
-#     @test collect(H_e*v_e) ≈ q_y.*v_e
-#     @test collect(H_s*v_s) ≈ q_x.*v_s
-#     @test collect(H_n*v_n) ≈ q_x.*v_n
-#
-#     @test collect(H_w'*v_w) == collect(H_w'*v_w)
-#     @test collect(H_e'*v_e) == collect(H_e'*v_e)
-#     @test collect(H_s'*v_s) == collect(H_s'*v_s)
-#     @test collect(H_n'*v_n) == collect(H_n'*v_n)
-# end
--- a/test/Grids/equidistant_grid_test.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -157,6 +157,9 @@
         ps = HyperBox((0,0),(2,1))
 
         @test equidistant_grid(ps, 3,4) == equidistant_grid((0,0), (2,1), 3,4)
+
+        @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
 
 
--- a/test/Grids/manifolds_test.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/test/Grids/manifolds_test.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -11,20 +11,31 @@
     @test ndims(unittetrahedron()) == 3
 end
 
+@testset "Interval" begin
+    @test Interval <: ParameterSpace{1}
+
+    @test Interval(0,1) isa Interval{Int}
+    @test Interval(0,1.) isa Interval{Float64}
+
+    @test unitinterval() isa Interval{Float64}
+    @test unitinterval() == Interval(0.,1.)
+    @test limits(unitinterval()) == (0.,1.)
+
+    @test unitinterval(Int) isa Interval{Int}
+    @test unitinterval(Int) == Interval(0,1)
+    @test limits(unitinterval(Int)) == (0,1)
+end
+
 @testset "HyperBox" begin
     @test HyperBox{<:Any, 2} <: ParameterSpace{2}
     @test HyperBox([1,1], [2,2]) isa HyperBox{Int, 2}
 
+    @test HyperBox([1,2], [1.,2.]) isa HyperBox{Float64,2}
+
     @test limits(HyperBox([1,2], [3,4])) == ([1,2], [3,4])
     @test limits(HyperBox([1,2], [3,4]), 1) == (1,3)
     @test limits(HyperBox([1,2], [3,4]), 2) == (2,4)
 
-    @test unitinterval() isa HyperBox{Float64,1}
-    @test limits(unitinterval()) == ([0], [1])
-
-    @test unitinterval(Int) isa HyperBox{Int,1}
-    @test limits(unitinterval(Int)) == ([0], [1])
-
     @test unitsquare() isa HyperBox{Float64,2}
     @test limits(unitsquare()) == ([0,0],[1,1])
 
@@ -40,6 +51,8 @@
     @test Simplex([1,2], [3,4]) isa Simplex{Int, 2}
     @test Simplex([1,2,3], [4,5,6],[1,1,1]) isa Simplex{Int, 3}
 
+    @test Simplex([1,2], [3.,4.]) isa Simplex{Float64, 2}
+
     @test verticies(Simplex([1,2], [3,4])) == ([1,2], [3,4])
 
     @test unittriangle() isa Simplex{Float64,2}
@@ -56,6 +69,7 @@
     @test c isa Chart{2}
     @test c([3,2]) == [6,4]
     @test parameterspace(c) == unitsquare()
+    @test ndims(c) == 2
 end
 
 @testset "Atlas" begin
--- a/test/Grids/mapped_grid_test.jl	Wed Sep 11 16:26:19 2024 +0200
+++ b/test/Grids/mapped_grid_test.jl	Mon Sep 16 10:33:47 2024 +0200
@@ -37,11 +37,11 @@
 
         @test mg isa Grid{SVector{2, Float64},2}
         @test jacobian(mg) isa Array{<:AbstractMatrix}
-        @test logicalgrid(mg) isa Grid
+        @test logical_grid(mg) isa Grid
 
         @test collect(mg) == x̄
         @test jacobian(mg) == J
-        @test logicalgrid(mg) == lg
+        @test logical_grid(mg) == lg
 
 
         x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
@@ -50,11 +50,11 @@
 
         @test mg isa Grid{SVector{3, Float64},2}
         @test jacobian(mg) isa Array{<:AbstractMatrix}
-        @test logicalgrid(mg) isa Grid
+        @test logical_grid(mg) isa Grid
 
         @test collect(mg) == x̄
         @test jacobian(mg) == J
-        @test logicalgrid(mg) == lg
+        @test logical_grid(mg) == lg
 
         sz1 = (10,11)
         sz2 = (10,12)
@@ -247,7 +247,7 @@
         ]
 
         function expected_bg(mg, bId, Jb)
-            lg = logicalgrid(mg)
+            lg = logical_grid(mg)
             return MappedGrid(
                 boundary_grid(lg, bId),
                 map(x̄, boundary_grid(lg, bId)),
@@ -279,34 +279,12 @@
     @test mg isa MappedGrid{SVector{2,Float64}, 2}
 
     lg = equidistant_grid((0,0), (1,1), 10, 11)
-    @test logicalgrid(mg) == lg
+    @test logical_grid(mg) == lg
     @test collect(mg) == map(x̄, lg)
 
     @test mapped_grid(lg, x̄, J) == mg
 end
 
-@testset "jacobian_determinant" begin
-    x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
-    J((ξ, η)) = @SMatrix[
-        η    ξ;
-        1   2η;
-    ]
-
-    g = mapped_grid(x̄, J, 10, 11)
-    J = map(logicalgrid(g)) do (ξ,η)
-        2η^2 - ξ
-    end
-    @test jacobian_determinant(g) ≈ J
-
-
-    lg = equidistant_grid((0,0), (1,1), 11, 21)
-    x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
-    J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg)
-    mg = MappedGrid(lg, x̄, J)
-
-    @test_broken jacobian(mg) isa AbstractArray{2,Float64}
-end
-
 @testset "metric_tensor" begin
     x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
     J((ξ, η)) = @SMatrix[
@@ -315,7 +293,7 @@
     ]
 
     g = mapped_grid(x̄, J, 10, 11)
-    G = map(logicalgrid(g)) do (ξ,η)
+    G = map(logical_grid(g)) do (ξ,η)
         @SMatrix[
             1+η^2   ξ*η+2η;
             ξ*η+2η  ξ^2 + 4η^2;
@@ -332,7 +310,7 @@
     ]
 
     g = mapped_grid(x̄, J, 10, 11)
-    G⁻¹ = map(logicalgrid(g)) do (ξ,η)
+    G⁻¹ = map(logical_grid(g)) do (ξ,η)
         @SMatrix[
             (1+η)^2  -ξ*(1+η);
             -ξ*(1+η) (1+ξ)^2+ξ^2;
@@ -384,7 +362,7 @@
     @test normal(g, CartesianBoundary{1,LowerBoundary}()) == fill(@SVector[-1,0], 11)
     @test normal(g, CartesianBoundary{1,UpperBoundary}()) == fill(@SVector[1,0], 11)
     @test normal(g, CartesianBoundary{2,LowerBoundary}()) == fill(@SVector[0,-1], 10)
-    @test normal(g, CartesianBoundary{2,UpperBoundary}()) ≈ map(boundary_grid(g,CartesianBoundary{2,UpperBoundary}())|>logicalgrid) do ξ̄
+    @test normal(g, CartesianBoundary{2,UpperBoundary}()) ≈ map(boundary_grid(g,CartesianBoundary{2,UpperBoundary}())|>logical_grid) do ξ̄
         α = 1-2ξ̄[1]
         @SVector[α,1]/√(α^2 + 1)
     end
@@ -393,28 +371,28 @@
 
     unit(v) = v/norm(v)
     @testset let bId = CartesianBoundary{1,LowerBoundary}()
-        lbg = boundary_grid(logicalgrid(g), bId)
+        lbg = boundary_grid(logical_grid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             -unit(@SVector[1/2,  η/3-1/6])
         end
     end
 
     @testset let bId = CartesianBoundary{1,UpperBoundary}()
-        lbg = boundary_grid(logicalgrid(g), bId)
+        lbg = boundary_grid(logical_grid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             unit(@SVector[7/2, 2η-1]/(5 + 3η + 2η^2))
         end
     end
 
     @testset let bId = CartesianBoundary{2,LowerBoundary}()
-        lbg = boundary_grid(logicalgrid(g), bId)
+        lbg = boundary_grid(logical_grid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ))
         end
     end
 
     @testset let bId = CartesianBoundary{2,UpperBoundary}()
-        lbg = boundary_grid(logicalgrid(g), bId)
+        lbg = boundary_grid(logical_grid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ))
         end
--- a/test/StaticDicts/StaticDicts_test.jl	Wed Sep 11 16:26:19 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,68 +0,0 @@
-using Test
-using Diffinitive.StaticDicts
-
-@testset "StaticDicts" begin
-
-@testset "StaticDict" begin
-    @testset "constructor" begin
-        @test (StaticDict{Int,Int,N} where N) <: AbstractDict
-
-        d = StaticDict(1=>2, 3=>4)
-        @test d isa StaticDict{Int,Int}
-        @test d[1] == 2
-        @test d[3] == 4
-
-        @test StaticDict((1=>2, 3=>4)) == d
-
-        @test StaticDict() isa StaticDict
-        @test StaticDict{Int,String}() isa StaticDict{Int,String,0}
-
-        @test StaticDict(1=>3, 2=>4.) isa StaticDict{Int,Real}
-        @test StaticDict(1. =>3, 2=>4) isa StaticDict{Real,Int}
-        @test StaticDict(1. =>3, 2=>4.) isa StaticDict{Real,Real}
-
-        @test_throws DomainError StaticDict(1=>3, 1=>3)
-    end
-
-    @testset "length" begin
-        @test length(StaticDict()) == 0
-        @test length(StaticDict(1=>1)) == 1
-        @test length(StaticDict(1=>1, 2=>2)) == 2
-    end
-
-    @testset "equality" begin
-        @test StaticDict(1=>1) == StaticDict(1=>1)
-        @test StaticDict(2=>1) != StaticDict(1=>1)
-        @test StaticDict(1=>2) != StaticDict(1=>1)
-
-        @test StaticDict(1=>1) === StaticDict(1=>1) #not true for a regular Dict
-        @test StaticDict(2=>1) !== StaticDict(1=>1)
-        @test StaticDict(1=>2) !== StaticDict(1=>1)
-    end
-
-    @testset "get" begin
-        d = StaticDict(1=>2, 3=>4)
-
-        @test get(d,1,6) == 2
-        @test get(d,3,6) == 4
-        @test get(d,5,6) == 6
-    end
-
-    @testset "iterate" begin
-        pairs = [1=>2, 3=>4, 5=>6]
-
-        d = StaticDict(pairs...)
-        @test collect(d) == pairs
-    end
-
-    @testset "merge" begin
-        @test merge(
-            StaticDict(1=>3, 2=> 4),
-            StaticDict(3=>5,4=>6)) == StaticDict(
-                1=>3, 2=>4, 3=>5, 4=>6
-            )
-        @test_throws DomainError merge(StaticDict(1=>3),StaticDict(1=>3))
-    end
-end
-
-end