--- a/Manifest.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/Manifest.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.4"
+julia_version = "1.10.5"
 manifest_format = "2.0"
 project_hash = "32fac879810099260f177c27318d3f26de4a00cc"
 
@@ -16,9 +16,9 @@
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "f54c23a5d304fb87110de62bace7777d59088c34"
+git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.15.0"
+version = "7.16.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -184,4 +184,4 @@
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.11.0+0"

--- a/Project.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/Project.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,7 +1,7 @@
-name = "Sbplib"
+name = "Diffinitive"
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-authors = ["Jonatan Werpers <jonatan@werpers.com>", "Vidar Stiernström <vidar.stiernstrom@it.uu.se>, and contributors"]
-version = "0.1.1"
+authors = ["Jonatan Werpers <jonatan@werpers.com>", "Vidar Stiernström <vidar.stiernstrom@gmail.com>, and contributors"]
+version = "0.1.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -12,10 +12,15 @@
 [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"
+Tokens = "040c2ec2-8d69-4aca-bf03-7d3a7092f2f6"
 
 [extensions]
-SbplibPlotsExt = "Plots"
-SbplibMakieExt = "Makie"
+DiffinitiveMakieExt = "Makie"
+DiffinitivePlotsExt = "Plots"
+DiffinitiveSparseArrayKitExt = ["SparseArrayKit", "Tokens"]
+DiffinitiveSparseArraysExt = ["SparseArrays", "Tokens"]
 
 [compat]
 julia = "1.5"

--- a/README.md	Wed Aug 28 10:35:08 2024 +0200
+++ b/README.md	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
-# Sbplib
+# Diffinitive
 
 ## Running tests
 To run all tests simply run
 ```
 (@v1.5) pkg> activate .
-(Sbplib) pkg> test
+(Diffinitive) pkg> test
 ```
 
 If you want to run tests from a specific file in `test/`, you can do
@@ -59,8 +59,8 @@
 
 ```julia
 using PkgBenchmark
-import Sbplib
-r = benchmarkpkg(Sbplib)
+import Diffinitive
+r = benchmarkpkg(Diffinitive)
 
 export_markdown(stdout, r)
 ```

--- a/TODO.md	Wed Aug 28 10:35:08 2024 +0200
+++ b/TODO.md	Wed Sep 11 15:41:58 2024 +0200
@@ -33,6 +33,7 @@
  - [ ] Kolla att vi har @inline på rätt ställen
  - [ ] Profilera
 
+ - [ ] Keep a lookout for allowing dependencies of package extensions (https://github.com/JuliaLang/Pkg.jl/issues/3641) This should be used to simplify the matrix extensions so that you don't have to load Tokens which is only used internally to the extension
 
 ### Grids
 

--- a/benchmark/Manifest.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/benchmark/Manifest.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.4"
+julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "25bba7b4a00465d5a2b00b589eb10e3301c31f2a"
+project_hash = "ecfc3e12aca5be17a874aba6134ff821abf61540"
 
 [[deps.AbstractTrees]]
 git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
@@ -25,9 +25,9 @@
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "f54c23a5d304fb87110de62bace7777d59088c34"
+git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.15.0"
+version = "7.16.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -99,6 +99,23 @@
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.Diffinitive]]
+deps = ["StaticArrays", "TOML", "TiledIteration"]
+path = ".."
+uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
+version = "0.1.1"
+
+    [deps.Diffinitive.extensions]
+    DiffinitiveMakieExt = "Makie"
+    DiffinitiveSparseArrayKitExt = ["SparseArrayKit", "Tokens"]
+    DiffinitiveSparseArraysExt = ["SparseArrays", "Tokens"]
+
+    [deps.Diffinitive.weakdeps]
+    Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+    SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    Tokens = "040c2ec2-8d69-4aca-bf03-7d3a7092f2f6"
+
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -273,18 +290,6 @@
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
 
-[[deps.Sbplib]]
-deps = ["StaticArrays", "TOML", "TiledIteration"]
-path = ".."
-uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-version = "0.1.1"
-
-    [deps.Sbplib.extensions]
-    SbplibMakieExt = "Makie"
-
-    [deps.Sbplib.weakdeps]
-    Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
@@ -391,7 +396,7 @@
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.11.0+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]

--- a/benchmark/Project.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/benchmark/Project.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Diffinitive = "5a373a26-915f-4769-bcab-bf03835de17b"
 Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
-PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
-Sbplib = "5a373a26-915f-4769-bcab-bf03835de17b"
+PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
\ No newline at end of file

--- a/benchmark/benchmark_laplace.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/benchmark/benchmark_laplace.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,7 +1,7 @@
-using Sbplib
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.RegionIndices
+using Diffinitive
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
 using BenchmarkTools
 
 # TODO: Move the below benchmarks into the benchmark suite

--- a/benchmark/benchmark_utils.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/benchmark/benchmark_utils.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -3,11 +3,11 @@
 import Mustache
 import Dates
 
-import Sbplib
+import Diffinitive
 
-const sbplib_root = splitpath(pathof(Sbplib))[1:end-2] |> joinpath
-const results_dir = mkpath(joinpath(sbplib_root, "benchmark/results"))
-const template_path = joinpath(sbplib_root, "benchmark/result.tmpl")
+const diffinitive_root = splitpath(pathof(Diffinitive))[1:end-2] |> joinpath
+const results_dir = mkpath(joinpath(diffinitive_root, "benchmark/results"))
+const template_path = joinpath(diffinitive_root, "benchmark/result.tmpl")
 
 """
     mainmain(;rev=nothing, target=nothing, baseline=nothing , kwargs...)
@@ -49,7 +49,7 @@
 `PkgBenchmark.BenchmarkResult`
 """
 function run_benchmark(;kwargs...)
-    r = PkgBenchmark.benchmarkpkg(Sbplib; kwargs...)
+    r = PkgBenchmark.benchmarkpkg(Diffinitive; kwargs...)
 
     rev = hg_rev() # Should be changed to hg_id() when the html can handle it.
 
@@ -158,17 +158,17 @@
 
 
 function hg_id()
-    cmd = Cmd(`hg id`, dir=sbplib_root)
+    cmd = Cmd(`hg id`, dir=diffinitive_root)
     return readchomp(addenv(cmd, "HGPLAIN"=>""))
 end
 
 function hg_rev()
-    cmd = Cmd(`hg id -i`, dir=sbplib_root)
+    cmd = Cmd(`hg id -i`, dir=diffinitive_root)
     return readchomp(addenv(cmd, "HGPLAIN"=>""))
 end
 
 function hg_update(rev)
-    cmd = Cmd(`hg update --check -r $rev`, dir=sbplib_root)
+    cmd = Cmd(`hg update --check -r $rev`, dir=diffinitive_root)
     run(addenv(cmd, "HGPLAIN"=>""))
 
     return nothing
@@ -182,9 +182,9 @@
 """
 function hg_commit(msg; secret=false)
     if secret
-        cmd = Cmd(`hg commit --verbose --secret --message $msg`, dir=sbplib_root)
+        cmd = Cmd(`hg commit --verbose --secret --message $msg`, dir=diffinitive_root)
     else
-        cmd = Cmd(`hg commit --verbose          --message $msg`, dir=sbplib_root)
+        cmd = Cmd(`hg commit --verbose          --message $msg`, dir=diffinitive_root)
     end
 
     out = readchomp(addenv(cmd, "HGPLAIN"=>""))
@@ -200,9 +200,9 @@
 """
 function hg_strip(rev; keep=false)
     if keep
-        cmd = Cmd(`hg --config extensions.strip= strip --keep -r $rev`, dir=sbplib_root)
+        cmd = Cmd(`hg --config extensions.strip= strip --keep -r $rev`, dir=diffinitive_root)
     else
-        cmd = Cmd(`hg --config extensions.strip= strip        -r $rev`, dir=sbplib_root)
+        cmd = Cmd(`hg --config extensions.strip= strip        -r $rev`, dir=diffinitive_root)
     end
 
     run(addenv(cmd, "HGPLAIN"=>""))
@@ -216,7 +216,7 @@
 Return true if the repositopry has uncommited changes.
 """
 function hg_is_dirty()
-    cmd = Cmd(`hg identify --id`, dir=sbplib_root)
+    cmd = Cmd(`hg identify --id`, dir=diffinitive_root)
     out = readchomp(addenv(cmd, "HGPLAIN"=>""))
 
     return endswith(out, "+")

--- a/benchmark/benchmarks.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/benchmark/benchmarks.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,9 @@
 using BenchmarkTools
 
-using Sbplib
-using Sbplib.Grids
-using Sbplib.SbpOperators
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
+using Diffinitive
+using Diffinitive.Grids
+using Diffinitive.SbpOperators
+using Diffinitive.LazyTensors
 
 using LinearAlgebra
 
@@ -168,20 +167,20 @@
 Dxx = second_derivative(g2, stencil_set, 1)
 Dyy = second_derivative(g2, stencil_set, 2)
 
-e₁ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,Lower}())
-e₁ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,Upper}())
-e₂ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,Lower}())
-e₂ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,Upper}())
+e₁ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,LowerBoundary}())
+e₁ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,UpperBoundary}())
+e₂ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,LowerBoundary}())
+e₂ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,UpperBoundary}())
 
-d₁ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{1,Lower}())
-d₁ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{1,Upper}())
-d₂ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{2,Lower}())
-d₂ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{2,Upper}())
+d₁ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{1,LowerBoundary}())
+d₁ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{1,UpperBoundary}())
+d₂ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{2,LowerBoundary}())
+d₂ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{2,UpperBoundary}())
 
-H₁ₗ = inner_product(boundary_grid(g2, CartesianBoundary{1,Lower}()), stencil_set)
-H₁ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{1,Upper}()), stencil_set)
-H₂ₗ = inner_product(boundary_grid(g2, CartesianBoundary{2,Lower}()), stencil_set)
-H₂ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{2,Upper}()), stencil_set)
+H₁ₗ = inner_product(boundary_grid(g2, CartesianBoundary{1,LowerBoundary}()), stencil_set)
+H₁ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{1,UpperBoundary}()), stencil_set)
+H₂ₗ = inner_product(boundary_grid(g2, CartesianBoundary{2,LowerBoundary}()), stencil_set)
+H₂ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{2,UpperBoundary}()), stencil_set)
 
 SUITE["boundary_terms"]["pre_composition"] = @benchmarkable $u2 .= $(H⁻¹∘e₁ₗ'∘H₁ₗ∘d₁ₗ)*$v2
 SUITE["boundary_terms"]["composition"]     = @benchmarkable $u2 .= ($H⁻¹∘$e₁ₗ'∘$H₁ₗ∘$d₁ₗ)*$v2

--- a/docs/Manifest.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/docs/Manifest.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.4"
+julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "4f0756199bb5f6739a5f4697152617efc4e0705c"
+project_hash = "c04450bb5c379e77d137cc05c4c0ab58eb1bfae9"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -30,9 +30,9 @@
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "f54c23a5d304fb87110de62bace7777d59088c34"
+git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.15.0"
+version = "7.16.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -94,6 +94,23 @@
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.Diffinitive]]
+deps = ["StaticArrays", "TOML", "TiledIteration"]
+path = ".."
+uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
+version = "0.1.1"
+
+    [deps.Diffinitive.extensions]
+    DiffinitiveMakieExt = "Makie"
+    DiffinitiveSparseArrayKitExt = ["SparseArrayKit", "Tokens"]
+    DiffinitiveSparseArraysExt = ["SparseArrays", "Tokens"]
+
+    [deps.Diffinitive.weakdeps]
+    Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+    SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    Tokens = "040c2ec2-8d69-4aca-bf03-7d3a7092f2f6"
+
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
 git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
@@ -102,9 +119,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 = "9d29b99b6b2b6bc2382a4c8dbec6eb694f389853"
+git-tree-sha1 = "5a1ee886566f2fa9318df1273d8b778b9d42712d"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.6.0"
+version = "1.7.0"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -149,9 +166,9 @@
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
+git-tree-sha1 = "f389674c99bfcde17dc57454011aa44d5a260a40"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.5.0"
+version = "1.6.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
@@ -246,9 +263,9 @@
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "a028ee3cb5641cccc4c24e90c36b0a4f7707bdf5"
+git-tree-sha1 = "1b35263570443fdd9e76c76b7062116e2f374ab8"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.0.14+0"
+version = "3.0.15+0"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -306,18 +323,6 @@
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
 
-[[deps.Sbplib]]
-deps = ["StaticArrays", "TOML", "TiledIteration"]
-path = ".."
-uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-version = "0.1.1"
-
-    [deps.Sbplib.extensions]
-    SbplibMakieExt = "Makie"
-
-    [deps.Sbplib.weakdeps]
-    Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
@@ -400,7 +405,7 @@
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.11.0+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]

--- a/docs/Project.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/docs/Project.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,3 +1,3 @@
 [deps]
+Diffinitive = "5a373a26-915f-4769-bcab-bf03835de17b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Sbplib = "5a373a26-915f-4769-bcab-bf03835de17b"

--- a/docs/make.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/docs/make.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,14 +1,14 @@
 using Documenter
-using Sbplib
+using Diffinitive
 
-using Sbplib.DiffOps
-using Sbplib.Grids
-using Sbplib.LazyTensors
-using Sbplib.RegionIndices
-using Sbplib.SbpOperators
-using Sbplib.StaticDicts
+using Diffinitive.DiffOps
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
+using Diffinitive.RegionIndices
+using Diffinitive.SbpOperators
+using Diffinitive.StaticDicts
 
-sitename = "Sbplib.jl"
+sitename = "Diffinitive.jl"
 
 remotes = nothing
 edit_link = nothing
@@ -31,6 +31,7 @@
     "Home" => "index.md",
     "operator_file_format.md",
     "grids_and_grid_functions.md",
+    "matrix_and_tensor_representations.md",
     "Submodules" => [
         "submodules/grids.md",
         "submodules/diff_ops.md",

--- a/docs/src/index.md	Wed Aug 28 10:35:08 2024 +0200
+++ b/docs/src/index.md	Wed Sep 11 15:41:58 2024 +0200
@@ -1,4 +1,4 @@
-# Sbplib.jl
+# Diffinitive.jl
 
 ```@contents
 Depth = 1

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/matrix_and_tensor_representations.md	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,5 @@
+# Matrix and tensor representations
+
+Sparse matrix and sparse tensor representations of lazy tensors can be constructed by loading [Tokens.jl](http://) and one of SparseArrays.jl or [SparseArrayKit.jl](http://). Through package extensions the following methods `sparse(::LazyTensor)` and `SparseArray(::LazyTensor)` are provided.
+
+<!-- TODO figure out how to add the docstrings here --/>

--- a/docs/src/operator_file_format.md	Wed Aug 28 10:35:08 2024 +0200
+++ b/docs/src/operator_file_format.md	Wed Sep 11 15:41:58 2024 +0200
@@ -1,6 +1,6 @@
 # Operator file format
 
-The intention is that Sbplib.jl should be a general and extensible framework
+The intention is that Diffinitive.jl should be a general and extensible framework
 for working with finite difference methods. It therefore includes a set of
 tools for storing and sharing operator definitions as well as a set of widely
 used operators.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/DiffinitiveMakieExt.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,88 @@
+module DiffinitiveMakieExt
+
+using Diffinitive.Grids
+using Makie
+using StaticArrays
+
+
+function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    ps = map(Tuple, g)[:]
+    values = gf[:]
+    faces = Vector{NTuple{3,Int}}()
+
+    n,m = size(g)
+    Li = LinearIndices((1:n, 1:m))
+    for i ∈ 1:n-1, j = 1:m-1
+
+        # Add point in the middle of the patch to preserve symmetries
+        push!(ps, Tuple((g[i,j] + g[i+1,j] + g[i+1,j+1] + g[i,j+1])/4))
+        push!(values, (gf[i,j] + gf[i+1,j] + gf[i+1,j+1] + gf[i,j+1])/4)
+
+        push!(faces, (Li[i,j],     Li[i+1,j],   length(ps)))
+        push!(faces, (Li[i+1,j],   Li[i+1,j+1], length(ps)))
+        push!(faces, (Li[i+1,j+1], Li[i,j+1],   length(ps)))
+        push!(faces, (Li[i,j+1],   Li[i,j],     length(ps)))
+    end
+
+    verticies = permutedims(reinterpret(reshape,eltype(eltype(ps)), ps))
+    faces = permutedims(reinterpret(reshape,Int, faces))
+
+    return verticies, faces, values
+end
+
+
+## Grids
+
+Makie.convert_arguments(::Type{<:Scatter}, g::Grid) = (reshape(map(Point,g),:),) # (map(Point,collect(g)[:]),)
+function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,2})
+    M = collect(g)
+
+    function cat_with_NaN(a,b)
+        vcat(a,[@SVector[NaN,NaN]],b)
+    end
+
+    xlines = reduce(cat_with_NaN, eachrow(M))
+    ylines = reduce(cat_with_NaN, eachcol(M))
+
+    return (cat_with_NaN(xlines,ylines),)
+end
+
+Makie.plot!(plot::Plot(Grid{<:Any,2})) = lines!(plot, plot.attributes, plot[1])
+
+
+## Grid functions
+
+### 1D
+function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
+    (collect(g), gf)
+end
+
+function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
+    (collect(g), gf)
+end
+
+Makie.plot!(plot::Plot(Grid{<:Any,1}, AbstractArray{<:Any,1})) = lines!(plot, plot.attributes, plot[1], plot[2])
+
+### 2D
+function Makie.convert_arguments(::Type{<:Surface}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    (getindex.(g,1), getindex.(g,2), gf)
+end
+
+function Makie.plot!(plot::Plot(Grid{<:Any,2},AbstractArray{<:Any, 2}))
+    r = @lift verticies_and_faces_and_values($(plot[1]), $(plot[2]))
+    v,f,c = (@lift $r[1]), (@lift $r[2]), (@lift $r[3])
+    mesh!(plot, plot.attributes, v, f;
+        color=c,
+        shading = NoShading,
+    )
+end
+# TBD: Can we define `mesh` instead of the above function and then forward plot! to that?
+
+function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    ps = map(g,gf) do (x,y), z
+        @SVector[x,y,z]
+    end
+    (reshape(ps,:),)
+end
+
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/DiffinitivePlotsExt.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,72 @@
+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
+
+# 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
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/DiffinitiveSparseArrayKitExt.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,22 @@
+module DiffinitiveSparseArrayKitExt
+
+using Diffinitive
+using Diffinitive.LazyTensors
+
+using SparseArrayKit
+using Tokens
+
+"""
+    SparseArray(t::LazyTensor)
+
+The sparse tensor representation of `t` with range dimensions to the left and
+domain dimensions to the right. If `L` is a `LazyTensor` with range and
+domain dimension 2 and `v` a 2-tensor, then `A = SparseArray(t)` is
+constructed so that `∑ₖ∑ₗA[i,j,k,l]*v[k,l] == L*v`.
+"""
+function SparseArrayKit.SparseArray(t::LazyTensor)
+    v = ArrayToken(:v, domain_size(t)...)
+    return Tokens._to_tensor(t*v, range_size(t), domain_size(t))
+end
+
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/DiffinitiveSparseArraysExt.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,25 @@
+module DiffinitiveSparseArraysExt
+
+using Diffinitive
+using Diffinitive.LazyTensors
+
+using SparseArrays
+using Tokens
+
+"""
+    sparse(t::LazyTensor)
+
+The sparse matrix representation of `t`.
+
+If `L` is a `LazyTensor` and `v` a tensor, then `A = sparse(L)` is constructed
+so that `A*reshape(v,:) == reshape(L*v,:)`.
+"""
+function SparseArrays.sparse(t::LazyTensor)
+    v = ArrayToken(:v, prod(domain_size(t)))
+
+    v̄ = reshape(v,domain_size(t)...)
+    tv = reshape(t*v̄, :)
+    return Tokens._to_matrix(tv, prod(range_size(t)), prod(domain_size(t)))
+end
+
+end

--- a/ext/SbplibMakieExt.jl	Wed Aug 28 10:35:08 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,88 +0,0 @@
-module SbplibMakieExt
-
-using Sbplib.Grids
-using Makie
-using StaticArrays
-
-
-function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
-    ps = map(Tuple, g)[:]
-    values = gf[:]
-    faces = Vector{NTuple{3,Int}}()
-
-    n,m = size(g)
-    Li = LinearIndices((1:n, 1:m))
-    for i ∈ 1:n-1, j = 1:m-1
-
-        # Add point in the middle of the patch to preserve symmetries
-        push!(ps, Tuple((g[i,j] + g[i+1,j] + g[i+1,j+1] + g[i,j+1])/4))
-        push!(values, (gf[i,j] + gf[i+1,j] + gf[i+1,j+1] + gf[i,j+1])/4)
-
-        push!(faces, (Li[i,j],     Li[i+1,j],   length(ps)))
-        push!(faces, (Li[i+1,j],   Li[i+1,j+1], length(ps)))
-        push!(faces, (Li[i+1,j+1], Li[i,j+1],   length(ps)))
-        push!(faces, (Li[i,j+1],   Li[i,j],     length(ps)))
-    end
-
-    verticies = permutedims(reinterpret(reshape,eltype(eltype(ps)), ps))
-    faces = permutedims(reinterpret(reshape,Int, faces))
-
-    return verticies, faces, values
-end
-
-
-## Grids
-
-Makie.convert_arguments(::Type{<:Scatter}, g::Grid) = (reshape(map(Point,g),:),) # (map(Point,collect(g)[:]),)
-function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,2})
-    M = collect(g)
-
-    function cat_with_NaN(a,b)
-        vcat(a,[@SVector[NaN,NaN]],b)
-    end
-
-    xlines = reduce(cat_with_NaN, eachrow(M))
-    ylines = reduce(cat_with_NaN, eachcol(M))
-
-    return (cat_with_NaN(xlines,ylines),)
-end
-
-Makie.plot!(plot::Plot(Grid{<:Any,2})) = lines!(plot, plot.attributes, plot[1])
-
-
-## Grid functions
-
-### 1D
-function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
-    (collect(g), gf)
-end
-
-function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
-    (collect(g), gf)
-end
-
-Makie.plot!(plot::Plot(Grid{<:Any,1}, AbstractArray{<:Any,1})) = lines!(plot, plot.attributes, plot[1], plot[2])
-
-### 2D
-function Makie.convert_arguments(::Type{<:Surface}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
-    (getindex.(g,1), getindex.(g,2), gf)
-end
-
-function Makie.plot!(plot::Plot(Grid{<:Any,2},AbstractArray{<:Any, 2}))
-    r = @lift verticies_and_faces_and_values($(plot[1]), $(plot[2]))
-    v,f,c = (@lift $r[1]), (@lift $r[2]), (@lift $r[3])
-    mesh!(plot, plot.attributes, v, f;
-        color=c,
-        shading = NoShading,
-    )
-end
-# TBD: Can we define `mesh` instead of the above function and then forward plot! to that?
-
-function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
-    ps = map(g,gf) do (x,y), z
-        @SVector[x,y,z]
-    end
-    (reshape(ps,:),)
-end
-
-end

--- a/ext/SbplibPlotsExt.jl	Wed Aug 28 10:35:08 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,72 +0,0 @@
-module SbplibPlotsExt
-
-using Sbplib.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
-
-# 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/sbp.jl	Wed Aug 28 10:35:08 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,9 +0,0 @@
-module sbp
-
-using Sbplib.Grids
-using Sbplib.RegionIndices
-using Sbplib.SbpOperators
-using Sbplib.DiffOps
-
-include("TimeStepper.jl")
-end  # module

--- a/sbpPlot.jl	Wed Aug 28 10:35:08 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,12 +0,0 @@
-include("sbp.jl")
-using Makie
-import .sbp.Grid
-function plotgridfunction(grid::sbp.Grid.EquidistantGrid, gridfunction::AbstractArray)
-    if sbp.Grid.dimension(grid) == 1
-        plot(sbp.Grid.pointsalongdim(grid,1), gridfunction)
-    elseif sbp.Grid.dimension(grid) == 2
-        scene = surface(sbp.Grid.pointsalongdim(grid,1),sbp.Grid.pointsalongdim(grid,2), gridfunction)
-    else
-        error(string("Plot not implemented for dimension ", string(dimension(grid))))
-    end
-end

--- a/src/DiffOps/DiffOps.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/DiffOps/DiffOps.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,9 +1,9 @@
 module DiffOps
 
-using Sbplib.RegionIndices
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.RegionIndices
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
 """
     DiffOp

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Diffinitive.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,16 @@
+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 Aug 28 10:35:08 2024 +0200
+++ b/src/Grids/Grids.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,7 @@
 # TODO: Double check that the interfaces for indexing and iterating are fully implemented and tested for all grids.
 module Grids
 
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
+using Diffinitive.LazyTensors
 using StaticArrays
 using LinearAlgebra
 
@@ -50,6 +49,8 @@
 export BoundaryIdentifier
 export TensorGridBoundary
 export CartesianBoundary
+export LowerBoundary
+export UpperBoundary
 
 export TensorGrid
 export ZeroDimGrid
@@ -69,8 +70,6 @@
 export metric_tensor
 export metric_tensor_inverse
 
-abstract type BoundaryIdentifier end
-
 include("manifolds.jl")
 include("grid.jl")
 include("tensor_grid.jl")

--- a/src/Grids/equidistant_grid.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -49,12 +49,30 @@
 
 min_spacing(g::EquidistantGrid) = spacing(g)
 
+"""
+    LowerBoundary <: BoundaryIdentifier
 
-boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
-boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
-boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
-boundary_indices(g::EquidistantGrid, id::Lower) = (1,)
-boundary_indices(g::EquidistantGrid, id::Upper) = (length(g),)
+Boundary identifier for the the lower (left) boundary of a one-dimensional grid.
+
+See also: [`BoundaryIdentifier`](@ref)
+"""
+struct LowerBoundary <: BoundaryIdentifier end
+
+"""
+    UpperBoundary <: BoundaryIdentifier
+
+Boundary identifier for the the upper (right)  boundary of a one-dimensional grid.
+
+See also: [`BoundaryIdentifier`](@ref)
+"""
+struct UpperBoundary <: BoundaryIdentifier end
+
+
+boundary_identifiers(::EquidistantGrid) = (LowerBoundary(), UpperBoundary())
+boundary_grid(g::EquidistantGrid, id::LowerBoundary) = ZeroDimGrid(g[begin])
+boundary_grid(g::EquidistantGrid, id::UpperBoundary) = ZeroDimGrid(g[end])
+boundary_indices(g::EquidistantGrid, id::LowerBoundary) = (firstindex(g),)
+boundary_indices(g::EquidistantGrid, id::UpperBoundary) = (lastindex(g),)
 
 """
     refine(g::EquidistantGrid, r::Int)

--- a/src/Grids/grid.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/Grids/grid.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -109,6 +109,13 @@
 function coarsen end
 
 """
+   BoundaryIdentifier
+
+An identifier for a boundary of a grid.
+"""
+abstract type BoundaryIdentifier end
+
+"""
     boundary_identifiers(g::Grid)
 
 Identifiers for all the boundaries of `g`.

--- a/src/Grids/mapped_grid.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,11 +1,58 @@
+"""
+    MappedGrid{T,D} <: Grid{T,D}
+
+A grid defined by a coordinate mapping from a logical grid to some physical
+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).
+"""
 struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D}
     logicalgrid::GT
     physicalcoordinates::CT
     jacobian::JT
+
+    """
+        MappedGrid(logicalgrid, 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))
+            throw(ArgumentError("Sizes must match"))
+        end
+
+        if size(first(jacobian)) != (length(first(physicalcoordinates)),D)
+            @show size(first(jacobian))
+            @show (length(first(physicalcoordinates)),D)
+            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)
+    end
 end
 
+function Base.:(==)(a::MappedGrid, b::MappedGrid)
+    same_logicalgrid = logicalgrid(a) == logicalgrid(b)
+    same_coordinates = collect(a) == collect(b)
+    same_jacobian = jacobian(a) == jacobian(b)
+
+    return same_logicalgrid && same_coordinates && same_jacobian
+end
+
+"""
+    logicalgrid(g::MappedGrid)
+
+The logical grid of `g`.
+"""
+logicalgrid(g::MappedGrid) = g.logicalgrid
+
+"""
+    jacobian(g::MappedGrid)
+
+The Jacobian matrix of `g` as a grid function.
+"""
 jacobian(g::MappedGrid) = g.jacobian
-logicalgrid(g::MappedGrid) = g.logicalgrid
 
 
 # Indexing interface
@@ -14,10 +61,8 @@
 
 Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d)
 Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d)
-# TODO: axes(...)?
 
 # Iteration interface
-
 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates)
 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state)
 
@@ -49,33 +94,65 @@
     )
 end
 
-# TBD: refine and coarsen could be implemented once we have a simple manifold implementation.
-# Before we do, we should consider the overhead of including such a field in the mapped grid struct.
+
+"""
+    mapped_grid(x, J, size::Vararg{Int})
 
-function mapped_grid(x, J, size...)
+A `MappedGrid` with a default logical grid on a unit hyper box. `x` and `J`
+are functions to be evaluated on the logical grid and `size` determines the
+size of the logical grid.
+"""
+function mapped_grid(x, J, size::Vararg{Int})
     D = length(size)
     lg = equidistant_grid(ntuple(i->0., D), ntuple(i->1., D), size...)
+    return mapped_grid(lg, x, J)
+end
+
+"""
+    mapped_grid(lg::Grid, x, J)
+
+A `MappedGrid` with logical grid `lg`. Physical coordinates and Jacobian are
+determined by the functions `x` and `J`.
+"""
+function mapped_grid(lg::Grid, x, J)
     return MappedGrid(
         lg,
         map(x,lg),
         map(J,lg),
     )
 end
-# TODO: Delete this function
 
+"""
+    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.
+"""
 function metric_tensor(g::MappedGrid)
     return map(jacobian(g)) do ∂x∂ξ
         ∂x∂ξ'*∂x∂ξ
     end
 end
 
+"""
+    metric_tensor_inverse(g::MappedGrid)
+
+The inverse of the metric tensor of `g` as a grid function.
+"""
 function metric_tensor_inverse(g::MappedGrid)
     return map(jacobian(g)) do ∂x∂ξ
         inv(∂x∂ξ'*∂x∂ξ)
@@ -119,7 +196,7 @@
 """
     normal(g::MappedGrid, boundary)
 
-The outward pointing normal as a grid function on the boundary
+The outward pointing normal as a grid function on the corresponding boundary grid.
 """
 function normal(g::MappedGrid, boundary)
     b_indices = boundary_indices(g, boundary)
@@ -132,11 +209,11 @@
 end
 
 function _boundary_sign(T, boundary)
-    if boundary_id(boundary) == Upper()
+    if boundary_id(boundary) == UpperBoundary()
         return one(T)
-    elseif boundary_id(boundary) == Lower()
+    elseif boundary_id(boundary) == LowerBoundary()
         return -one(T)
     else
-        throw(ArgumentError("The boundary identifier must be either `Lower()` or `Upper()`"))
+        throw(ArgumentError("The boundary identifier must be either `LowerBoundary()` or `UpperBoundary()`"))
     end
 end

--- a/src/Grids/tensor_grid.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/Grids/tensor_grid.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,3 @@
-# TODO: Check this file and other grids for duplicate implementation of general methods implemented for Grid
-
 """
     TensorGrid{T,D} <: Grid{T,D}
 
@@ -119,7 +117,7 @@
 end
 
 """
-   grid_and_local_dim_index(nds, d)
+    grid_and_local_dim_index(nds, d)
 
 Given a tuple of number of dimensions `nds`, and a global dimension index `d`,
 calculate which grid index, and local dimension, `d` corresponds to.

--- a/src/SbpOperators/SbpOperators.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -40,9 +40,9 @@
 export sat_tensors
 
 # Using
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
-using Sbplib.Grids
+using Diffinitive.RegionIndices
+using Diffinitive.LazyTensors
+using Diffinitive.Grids
 
 # Includes
 include("stencil.jl")

--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,27 +1,27 @@
 """
-    BoundaryOperator{T,R,N} <: LazyTensor{T,0,1}
+    BoundaryOperator{T,B,N} <: LazyTensor{T,0,1}
 
 Implements the boundary operator `op` for 1D as a `LazyTensor`
 
 `op` is the restriction of a grid function to the boundary using some closure
-`Stencil{T,N}`. The boundary to restrict to is determined by `R`. `op'` is the
+`Stencil{T,N}`. The boundary to restrict to is determined by `B`. `op'` is the
 prolongation of a zero dimensional array to the whole grid using the same
 closure stencil.
 """
-struct BoundaryOperator{T,R<:Region,N} <: LazyTensor{T,0,1}
+struct BoundaryOperator{T,B<:BoundaryIdentifier,N} <: LazyTensor{T,0,1}
     stencil::Stencil{T,N}
     size::Int
 end
 
 """
-    BoundaryOperator(grid::EquidistantGrid, closure_stencil, region)
+    BoundaryOperator(grid::EquidistantGrid, closure_stencil, boundary)
 
 Constructs the BoundaryOperator with stencil `closure_stencil` for a
 `EquidistantGrid` `grid`, restricting to to the boundary specified by
-`region`.
+`boundary`.
 """
-function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, region::Region) where {T,N}
-    return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1])
+function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, boundary::BoundaryIdentifier) where {T,N}
+    return BoundaryOperator{T,typeof(boundary),N}(closure_stencil,size(grid)[1])
 end
 
 """
@@ -29,24 +29,24 @@
 
 The size of the closure stencil.
 """
-closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N
+closure_size(::BoundaryOperator{T,B,N}) where {T,B,N} = N
 
 LazyTensors.range_size(op::BoundaryOperator) = ()
 LazyTensors.domain_size(op::BoundaryOperator) = (op.size,)
 
-function LazyTensors.apply(op::BoundaryOperator{<:Any,Lower}, v::AbstractVector)
+function LazyTensors.apply(op::BoundaryOperator{<:Any,LowerBoundary}, v::AbstractVector)
     apply_stencil(op.stencil,v,1)
 end
 
-function LazyTensors.apply(op::BoundaryOperator{<:Any,Upper}, v::AbstractVector)
+function LazyTensors.apply(op::BoundaryOperator{<:Any,UpperBoundary}, v::AbstractVector)
     apply_stencil_backwards(op.stencil,v,op.size)
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Lower}, v::AbstractArray{<:Any,0}, i::Index{Lower})
+function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,LowerBoundary}, v::AbstractArray{<:Any,0}, i::Index{Lower})
     return op.stencil[Int(i)-1]*v[]
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Upper}, v::AbstractArray{<:Any,0}, i::Index{Upper})
+function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,UpperBoundary}, v::AbstractArray{<:Any,0}, i::Index{Upper})
     return op.stencil[op.size[1] - Int(i)]*v[]
 end
 

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -60,7 +60,7 @@
 condition. `H_tuning` and `R_tuning` are used to specify the strength of the
 penalty.
 
-See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
+See also: [`sat`](@ref), [`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
 """
 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
     id = boundary(bc)
@@ -69,7 +69,7 @@
     Hᵧ = inner_product(boundary_grid(g, id), set)
     e = boundary_restriction(g, set, id)
     d = normal_derivative(g, set, id)
-    B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning)
+    B = positivity_decomposition(Δ, g, boundary(bc); H_tuning, R_tuning)
     penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
     return penalty_tensor, e
 end
@@ -94,7 +94,7 @@
 end
 
 """
-    positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+    positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
 
 Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive
 definite with respect to the boundary quadrature. Here `d` is the normal
@@ -102,28 +102,26 @@
 to form a symmetric and energy stable penalty for a Dirichlet condition. The
 parameters `H_tuning` and `R_tuning` are used to specify the strength of the
 penalty and must be greater than 1. For details we refer to
-https://doi.org/10.1016/j.jcp.2020.109294
+<https://doi.org/10.1016/j.jcp.2020.109294>
 """
-function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+function positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
     @assert(H_tuning ≥ 1.)
     @assert(R_tuning ≥ 1.)
-    Nτ_H, τ_R = positivity_limits(Δ,g,bc)
+    Nτ_H, τ_R = positivity_limits(Δ,g,b)
     return H_tuning*Nτ_H + R_tuning*τ_R
 end
 
-# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
-# change bc::BoundaryCondition to id::BoundaryIdentifier
-function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
+function positivity_limits(Δ::Laplace, g::EquidistantGrid, b::BoundaryIdentifier)
     h = spacing(g)
     θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
     θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
 
-    τ_H = 1/(h*θ_H)
-    τ_R = 1/(h*θ_R)
+    τ_H = one(eltype(Δ))/(h*θ_H)
+    τ_R = one(eltype(Δ))/(h*θ_R)
     return τ_H, τ_R
 end
 
-function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition)
-    τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc)
+function positivity_limits(Δ::Laplace, g::TensorGrid, b::BoundaryIdentifier)
+    τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(b)], b)
     return τ_H*ndims(g), τ_R
 end

--- a/src/Sbplib.jl	Wed Aug 28 10:35:08 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,16 +0,0 @@
-module Sbplib
-
-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/test/DiffOps/DiffOps_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/DiffOps/DiffOps_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,9 +1,9 @@
 using Test
-using Sbplib.DiffOps
-using Sbplib.Grids
-using Sbplib.SbpOperators
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
+using Diffinitive.DiffOps
+using Diffinitive.Grids
+using Diffinitive.SbpOperators
+using Diffinitive.RegionIndices
+using Diffinitive.LazyTensors
 
 #
 # @testset "BoundaryValue" begin

--- a/test/Grids/equidistant_grid_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,7 +1,6 @@
-using Sbplib.Grids
+using Diffinitive.Grids
 using Test
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
+using Diffinitive.LazyTensors
 using StaticArrays
 
 
@@ -64,24 +63,24 @@
 
     @testset "boundary_identifiers" begin
         g = EquidistantGrid(0:0.1:10)
-        @test boundary_identifiers(g) == (Lower(), Upper())
+        @test boundary_identifiers(g) == (LowerBoundary(), UpperBoundary())
         @inferred boundary_identifiers(g)
     end
 
     @testset "boundary_grid" begin
         g = EquidistantGrid(0:0.1:1)
-        @test boundary_grid(g, Lower()) == ZeroDimGrid(0.0)
-        @test boundary_grid(g, Upper()) == ZeroDimGrid(1.0)
+        @test boundary_grid(g, LowerBoundary()) == ZeroDimGrid(0.0)
+        @test boundary_grid(g, UpperBoundary()) == ZeroDimGrid(1.0)
     end
 
     @testset "boundary_indices" begin
         g = EquidistantGrid(0:0.1:1)
-        @test boundary_indices(g, Lower()) == (1,)
-        @test boundary_indices(g, Upper()) == (11,)
+        @test boundary_indices(g, LowerBoundary()) == (1,)
+        @test boundary_indices(g, UpperBoundary()) == (11,)
 
         g = EquidistantGrid(2:0.1:10)
-        @test boundary_indices(g, Lower()) == (1,)
-        @test boundary_indices(g, Upper()) == (81,)
+        @test boundary_indices(g, LowerBoundary()) == (1,)
+        @test boundary_indices(g, UpperBoundary()) == (81,)
 
     end
 

--- a/test/Grids/grid_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Grids/grid_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,6 +1,6 @@
 using Test
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 using StaticArrays
 
 @testset "Grid" begin

--- a/test/Grids/manifolds_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Grids/manifolds_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 using Test
 
-using Sbplib.Grids
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
+using Diffinitive.LazyTensors
 
 # using StaticArrays
 

--- a/test/Grids/mapped_grid_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Grids/mapped_grid_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
-using Sbplib.Grids
-using Sbplib.RegionIndices
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
 using Test
 using StaticArrays
 using LinearAlgebra
@@ -28,33 +28,82 @@
 end
 
 @testset "MappedGrid" begin
-    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)
+    @testset "Constructor" begin
+        lg = equidistant_grid((0,0), (1,1), 11, 21)
+
+        x̄ = map(ξ̄ -> 2ξ̄, lg)
+        J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
+        mg = MappedGrid(lg, x̄, J)
 
-    # TODO: Test constructor for different dims of range and domain for the coordinates
-    # TODO: Test constructor with different type than TensorGrid. a dummy type?
+        @test mg isa Grid{SVector{2, Float64},2}
+        @test jacobian(mg) isa Array{<:AbstractMatrix}
+        @test logicalgrid(mg) isa Grid
 
-    @test_broken false # @test_throws ArgumentError("Sizes must match") MappedGrid(lg, map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11))
+        @test collect(mg) == x̄
+        @test jacobian(mg) == J
+        @test logicalgrid(mg) == lg
 
 
-    @test mg isa Grid{SVector{2, Float64},2}
+        x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
+        J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg)
+        mg = MappedGrid(lg, x̄, J)
+
+        @test mg isa Grid{SVector{3, Float64},2}
+        @test jacobian(mg) isa Array{<:AbstractMatrix}
+        @test logicalgrid(mg) isa Grid
+
+        @test collect(mg) == x̄
+        @test jacobian(mg) == J
+        @test logicalgrid(mg) == lg
+
+        sz1 = (10,11)
+        sz2 = (10,12)
+        @test_throws ArgumentError("Sizes must match") MappedGrid(
+            equidistant_grid((0,0), (1,1), sz2...),
+            rand(SVector{2},sz1...),
+            rand(SMatrix{2,2},sz1...),
+        )
 
-    @test jacobian(mg) isa Array{<:AbstractMatrix}
-    @test logicalgrid(mg) isa Grid
+        @test_throws ArgumentError("Sizes must match") MappedGrid(
+            equidistant_grid((0,0), (1,1), sz1...),
+            rand(SVector{2},sz2...),
+            rand(SMatrix{2,2},sz1...),
+        )
+
+        @test_throws ArgumentError("Sizes must match") MappedGrid(
+            equidistant_grid((0,0), (1,1), sz1...),
+            rand(SVector{2},sz1...),
+            rand(SMatrix{2,2},sz2...),
+        )
+
+        err_str = "The size of the jacobian must match the dimensions of the grid and coordinates"
+        @test_throws ArgumentError(err_str) MappedGrid(
+            equidistant_grid((0,0), (1,1), 10, 11),
+            rand(SVector{3}, 10, 11),
+            rand(SMatrix{3,4}, 10, 11),
+        )
+
+        @test_throws ArgumentError(err_str) MappedGrid(
+            equidistant_grid((0,0), (1,1), 10, 11),
+            rand(SVector{3}, 10, 11),
+            rand(SMatrix{4,2}, 10, 11),
+        )
+    end
 
     @testset "Indexing Interface" begin
+        lg = equidistant_grid((0,0), (1,1), 11, 21)
+        x̄ = map(ξ̄ -> 2ξ̄, lg)
+        J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
         mg = MappedGrid(lg, x̄, J)
         @test mg[1,1] == [0.0, 0.0]
-        @test mg[4,2] == [0.6, 0.2]
-        @test mg[6,10] == [1., 1.8]
+        @test mg[4,2] == [0.6, 0.1]
+        @test mg[6,10] == [1., 0.9]
 
         @test mg[begin, begin] == [0.0, 0.0]
         @test mg[end,end] == [2.0, 2.0]
         @test mg[begin,end] == [0., 2.]
 
-        @test eachindex(mg) == CartesianIndices((11,11))
+        @test axes(mg) == (1:11, 1:21)
 
         @testset "cartesian indexing" begin
             cases = [
@@ -71,7 +120,7 @@
         end
 
         @testset "eachindex" begin
-            @test eachindex(mg) == CartesianIndices((11,11))
+            @test eachindex(mg) == CartesianIndices((11,21))
         end
 
         @testset "firstindex" begin
@@ -81,15 +130,21 @@
 
         @testset "lastindex" begin
             @test lastindex(mg, 1) == 11
-            @test lastindex(mg, 2) == 11
+            @test lastindex(mg, 2) == 21
         end
     end
-    # TODO: Test with different types of logical grids
 
     @testset "Iterator interface" begin
+        lg = equidistant_grid((0,0), (1,1), 11, 21)
+        x̄ = map(ξ̄ -> 2ξ̄, lg)
+        J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
+
+        mg = MappedGrid(lg, x̄, J)
+
+        lg2 = equidistant_grid((0,0), (1,1), 15, 11)
         sg = MappedGrid(
             equidistant_grid((0,0), (1,1), 15, 11),
-            map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11)
+            map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg2), rand(SMatrix{3,2,Float64},15,11)
         )
 
         @test eltype(mg) == SVector{2,Float64}
@@ -98,13 +153,13 @@
         @test eltype(typeof(mg)) == SVector{2,Float64}
         @test eltype(typeof(sg)) == SVector{3,Float64}
 
-        @test size(mg) == (11,11)
+        @test size(mg) == (11,21)
         @test size(sg) == (15,11)
 
-        @test size(mg,2) == 11
+        @test size(mg,2) == 21
         @test size(sg,2) == 11
 
-        @test length(mg) == 121
+        @test length(mg) == 231
         @test length(sg) == 165
 
         @test Base.IteratorSize(mg) == Base.HasShape{2}()
@@ -127,17 +182,56 @@
     end
 
     @testset "Base" begin
+        lg = equidistant_grid((0,0), (1,1), 11, 21)
+        x̄ = map(ξ̄ -> 2ξ̄, lg)
+        J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
+        mg = MappedGrid(lg, x̄, J)
+
         @test ndims(mg) == 2
     end
 
+    @testset "==" begin
+        sz = (15,11)
+        lg = equidistant_grid((0,0), (1,1), sz...)
+        x = rand(SVector{3,Float64}, sz...)
+        J = rand(SMatrix{3,2,Float64}, sz...)
+
+        sg = MappedGrid(lg, x, J)
+
+        sg1 = MappedGrid(equidistant_grid((0,0), (1,1), sz...), copy(x), copy(J))
+
+        sz2 = (15,12)
+        lg2 = equidistant_grid((0,0), (1,1), sz2...)
+        x2 = rand(SVector{3,Float64}, sz2...)
+        J2 = rand(SMatrix{3,2,Float64}, sz2...)
+        sg2 = MappedGrid(lg2, x2, J2)
+
+        sg3 = MappedGrid(lg, rand(SVector{3,Float64}, sz...), J)
+        sg4 = MappedGrid(lg, x, rand(SMatrix{3,2,Float64}, sz...))
+
+        @test sg == sg1
+        @test sg != sg2 # Different size
+        @test sg != sg3 # Different coordinates
+        @test sg != sg4 # Different jacobian
+    end
+
     @testset "boundary_identifiers" begin
+        lg = equidistant_grid((0,0), (1,1), 11, 15)
+        x̄ = map(ξ̄ -> 2ξ̄, lg)
+        J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
+        mg = MappedGrid(lg, x̄, J)
         @test boundary_identifiers(mg) == boundary_identifiers(lg)
     end
 
     @testset "boundary_indices" begin
-        @test boundary_indices(mg, CartesianBoundary{1,Lower}()) == boundary_indices(lg,CartesianBoundary{1,Lower}())
-        @test boundary_indices(mg, CartesianBoundary{2,Lower}()) == boundary_indices(lg,CartesianBoundary{2,Lower}())
-        @test boundary_indices(mg, CartesianBoundary{1,Upper}()) == boundary_indices(lg,CartesianBoundary{1,Upper}())
+        lg = equidistant_grid((0,0), (1,1), 11, 15)
+        x̄ = map(ξ̄ -> 2ξ̄, lg)
+        J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
+        mg = MappedGrid(lg, x̄, J)
+
+        @test boundary_indices(mg, CartesianBoundary{1,LowerBoundary}()) == boundary_indices(lg,CartesianBoundary{1,LowerBoundary}())
+        @test boundary_indices(mg, CartesianBoundary{2,LowerBoundary}()) == boundary_indices(lg,CartesianBoundary{2,LowerBoundary}())
+        @test boundary_indices(mg, CartesianBoundary{1,UpperBoundary}()) == boundary_indices(lg,CartesianBoundary{1,UpperBoundary}())
     end
 
     @testset "boundary_grid" begin
@@ -152,28 +246,30 @@
             1+ξ*(ξ-1);
         ]
 
-        function test_boundary_grid(mg, bId, Jb)
-            bg = boundary_grid(mg, bId)
-
+        function expected_bg(mg, bId, Jb)
             lg = logicalgrid(mg)
-            expected_bg = MappedGrid(
+            return MappedGrid(
                 boundary_grid(lg, bId),
                 map(x̄, boundary_grid(lg, bId)),
                 map(Jb, boundary_grid(lg, bId)),
             )
+        end
 
-            @testset let bId=bId, bg=bg, expected_bg=expected_bg
-                @test collect(bg) == collect(expected_bg)
-                @test logicalgrid(bg) == logicalgrid(expected_bg)
-                @test jacobian(bg) == jacobian(expected_bg)
-                # TODO: Implement equality of a curvilinear grid and simlify the above
-            end
+        let bid = TensorGridBoundary{1, LowerBoundary}()
+            @test boundary_grid(mg, bid) == expected_bg(mg, bid, J2)
         end
 
-        @testset test_boundary_grid(mg, TensorGridBoundary{1, Lower}(), J2)
-        @testset test_boundary_grid(mg, TensorGridBoundary{1, Upper}(), J2)
-        @testset test_boundary_grid(mg, TensorGridBoundary{2, Lower}(), J1)
-        @testset test_boundary_grid(mg, TensorGridBoundary{2, Upper}(), J1)
+        let bid = TensorGridBoundary{1, UpperBoundary}()
+            @test boundary_grid(mg, bid) == expected_bg(mg, bid, J2)
+        end
+
+        let bid = TensorGridBoundary{2, LowerBoundary}()
+            @test boundary_grid(mg, bid) == expected_bg(mg, bid, J1)
+        end
+
+        let bid = TensorGridBoundary{2, UpperBoundary}()
+            @test boundary_grid(mg, bid) == expected_bg(mg, bid, J1)
+        end
     end
 end
 
@@ -185,18 +281,66 @@
     lg = equidistant_grid((0,0), (1,1), 10, 11)
     @test logicalgrid(mg) == lg
     @test collect(mg) == map(x̄, lg)
+
+    @test mapped_grid(lg, x̄, J) == mg
 end
 
 @testset "jacobian_determinant" begin
-    @test_broken false
+    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
-    @test_broken false
+    x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
+    J((ξ, η)) = @SMatrix[
+        η    ξ;
+        1   2η;
+    ]
+
+    g = mapped_grid(x̄, J, 10, 11)
+    G = map(logicalgrid(g)) do (ξ,η)
+        @SMatrix[
+            1+η^2   ξ*η+2η;
+            ξ*η+2η  ξ^2 + 4η^2;
+        ]
+    end
+    @test metric_tensor(g) ≈ G
 end
 
 @testset "metric_tensor_inverse" begin
-    @test_broken false
+    x̄((ξ, η)) = @SVector[ξ + ξ^2/2, η + η^2 + ξ^2/2]
+    J((ξ, η)) = @SMatrix[
+        1+ξ   0;
+        ξ    1+η;
+    ]
+
+    g = mapped_grid(x̄, J, 10, 11)
+    G⁻¹ = map(logicalgrid(g)) do (ξ,η)
+        @SMatrix[
+            (1+η)^2  -ξ*(1+η);
+            -ξ*(1+η) (1+ξ)^2+ξ^2;
+        ]/(((1+ξ)^2+ξ^2)*(1+η)^2 - ξ^2*(1+η)^2)
+
+    end
+
+    @test metric_tensor_inverse(g) ≈ G⁻¹
 end
 
 @testset "min_spacing" begin
@@ -237,10 +381,10 @@
 @testset "normal" begin
     g = mapped_grid(_partially_curved_mapping()...,10, 11)
 
-    @test normal(g, CartesianBoundary{1,Lower}()) == fill(@SVector[-1,0], 11)
-    @test normal(g, CartesianBoundary{1,Upper}()) == fill(@SVector[1,0], 11)
-    @test normal(g, CartesianBoundary{2,Lower}()) == fill(@SVector[0,-1], 10)
-    @test normal(g, CartesianBoundary{2,Upper}()) ≈ map(boundary_grid(g,CartesianBoundary{2,Upper}())|>logicalgrid) do ξ̄
+    @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 ξ̄
         α = 1-2ξ̄[1]
         @SVector[α,1]/√(α^2 + 1)
     end
@@ -248,28 +392,28 @@
     g = mapped_grid(_fully_curved_mapping()...,5,4)
 
     unit(v) = v/norm(v)
-    @testset let bId = CartesianBoundary{1,Lower}()
+    @testset let bId = CartesianBoundary{1,LowerBoundary}()
         lbg = boundary_grid(logicalgrid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             -unit(@SVector[1/2,  η/3-1/6])
         end
     end
 
-    @testset let bId = CartesianBoundary{1,Upper}()
+    @testset let bId = CartesianBoundary{1,UpperBoundary}()
         lbg = boundary_grid(logicalgrid(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,Lower}()
+    @testset let bId = CartesianBoundary{2,LowerBoundary}()
         lbg = boundary_grid(logicalgrid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ))
         end
     end
 
-    @testset let bId = CartesianBoundary{2,Upper}()
+    @testset let bId = CartesianBoundary{2,UpperBoundary}()
         lbg = boundary_grid(logicalgrid(g), bId)
         @test normal(g, bId) ≈ map(lbg) do (ξ, η)
             unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ))

--- a/test/Grids/tensor_grid_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Grids/tensor_grid_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,7 +1,6 @@
 using Test
-using Sbplib.Grids
+using Diffinitive.Grids
 using StaticArrays
-using Sbplib.RegionIndices
 
 @testset "TensorGrid" begin
     g₁ = EquidistantGrid(range(0,1,length=11))
@@ -170,13 +169,13 @@
     end
 
     @testset "boundary_identifiers" begin
-        @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (Lower,Upper,Lower,Upper))
-        @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,Lower}(),TensorGridBoundary{1,Upper}())
+        @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (LowerBoundary,UpperBoundary,LowerBoundary,UpperBoundary))
+        @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,LowerBoundary}(),TensorGridBoundary{1,UpperBoundary}())
     end
 
     @testset "boundary_grid" begin
-        @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂)
-        @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄)
+        @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, UpperBoundary}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂)
+        @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, UpperBoundary}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄)
     end
 
     @testset "boundary_indices" begin
@@ -184,14 +183,14 @@
         g₂ = EquidistantGrid(range(2,3,length=6))
         g₄ = ZeroDimGrid(@SVector[1,2])
 
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Lower}()) == (1,:)
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == (11,:)
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Lower}()) == (:,1)
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Upper}()) == (:,6)
-        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Lower}()) == (1,)
-        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == (11,)
-        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Lower}()) == (1,)
-        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Upper}()) == (11,)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, LowerBoundary}()) == (1,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, UpperBoundary}()) == (11,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, LowerBoundary}()) == (:,1)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, UpperBoundary}()) == (:,6)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, LowerBoundary}()) == (1,)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, UpperBoundary}()) == (11,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, LowerBoundary}()) == (1,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, UpperBoundary}()) == (11,)
     end
 end
 

--- a/test/Grids/zero_dim_grid_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Grids/zero_dim_grid_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
 using Test
-using Sbplib.Grids
+using Diffinitive.Grids
 using StaticArrays
 
 @testset "ZeroDimGrid" begin

--- a/test/LazyTensors/lazy_array_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/LazyTensors/lazy_array_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,6 +1,6 @@
 using Test
-using Sbplib.LazyTensors
-using Sbplib.RegionIndices
+using Diffinitive.LazyTensors
+using Diffinitive.RegionIndices
 
 
 @testset "LazyArray" begin

--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,6 +1,6 @@
 using Test
-using Sbplib.LazyTensors
-using Sbplib.RegionIndices
+using Diffinitive.LazyTensors
+using Diffinitive.RegionIndices
 
 using Tullio
 

--- a/test/LazyTensors/lazy_tensor_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/LazyTensors/lazy_tensor_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
 using Test
-using Sbplib.LazyTensors
+using Diffinitive.LazyTensors
 
 @testset "Generic Mapping methods" begin
     struct DummyMapping{T,R,D} <: LazyTensor{T,R,D} end

--- a/test/LazyTensors/tensor_types_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/LazyTensors/tensor_types_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
 using Test
-using Sbplib.LazyTensors
+using Diffinitive.LazyTensors
 using BenchmarkTools
 
 @testset "IdentityTensor" begin

--- a/test/LazyTensors/tuple_manipulation_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/LazyTensors/tuple_manipulation_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
 using Test
-using Sbplib.LazyTensors
+using Diffinitive.LazyTensors
 
 @testset "split_index" begin
     @test LazyTensors.split_index(2,1,2,2, 1,2,3,4,5,6) == ((1,2,:,5,6),(3,4))

--- a/test/Manifest.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Manifest.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.4"
+julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "f2b0634c12bbed93a17efc88d466604d5a07c465"
+project_hash = "9dddd5385164ee197d1b3f22302bc95701c1f5e5"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -74,9 +74,9 @@
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
+git-tree-sha1 = "f389674c99bfcde17dc57454011aa44d5a260a40"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.5.0"
+version = "1.6.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
@@ -176,6 +176,12 @@
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
+[[deps.PackageExtensionCompat]]
+git-tree-sha1 = "fb28e33b8a95c4cee25ce296c817d89cc2e53518"
+uuid = "65ce6f38-6b18-4e1d-a461-8949797d7930"
+version = "1.0.2"
+weakdeps = ["Requires", "TOML"]
+
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
 git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821"
@@ -231,6 +237,20 @@
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[deps.SparseArrayKit]]
+deps = ["LinearAlgebra", "PackageExtensionCompat", "TupleTools", "VectorInterface"]
+git-tree-sha1 = "e06f75c460c16b5b08bf93324b7db2cf9c1fe831"
+uuid = "a9a3c162-d163-4c15-8926-b8794fbefed2"
+version = "0.3.1"
+
+    [deps.SparseArrayKit.extensions]
+    SparseArrayKitSparseArrays = "SparseArrays"
+    SparseArrayKitTensorOperations = "TensorOperations"
+
+    [deps.SparseArrayKit.weakdeps]
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
+
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -297,6 +317,12 @@
 uuid = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 version = "3.0.0"
 
+[[deps.Tokens]]
+deps = ["SparseArrayKit", "SparseArrays"]
+git-tree-sha1 = "c4f40125383ce3bfcfcd49a1b206080b7afd9a34"
+uuid = "040c2ec2-8d69-4aca-bf03-7d3a7092f2f6"
+version = "0.1.1"
+
 [[deps.Tullio]]
 deps = ["DiffRules", "LinearAlgebra", "Requires"]
 git-tree-sha1 = "6d476962ba4e435d7f4101a403b1d3d72afe72f3"
@@ -315,6 +341,11 @@
     FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
+[[deps.TupleTools]]
+git-tree-sha1 = "41d61b1c545b06279871ef1a4b5fcb2cac2191cd"
+uuid = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
+version = "1.5.0"
+
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -322,6 +353,12 @@
 [[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
+[[deps.VectorInterface]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "7aff7d62bffad9bba9928eb6ab55226b32a351eb"
+uuid = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
+version = "0.4.6"
+
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
@@ -330,7 +367,7 @@
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.11.0+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]

--- a/test/Project.toml	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/Project.toml	Wed Sep 11 15:41:58 2024 +0200
@@ -2,8 +2,11 @@
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
+Tokens = "040c2ec2-8d69-4aca-bf03-7d3a7092f2f6"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

--- a/test/RegionIndices/RegionIndices_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/RegionIndices/RegionIndices_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,3 +1,3 @@
-using Sbplib.RegionIndices
+using Diffinitive.RegionIndices
 using Test
 

--- a/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 using Test
 
-using Sbplib.Grids
-using Sbplib.RegionIndices
-using Sbplib.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
+using Diffinitive.SbpOperators
 
 @testset "BoundaryCondition" begin
     grid_1d = equidistant_grid(0.0, 1.0, 11)
@@ -15,8 +15,8 @@
     g = 3.14
     f(x,y,z) = x^2+y^2+z^2
     @testset "Constructors" begin
-        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower}
-        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,Lower}}
+        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,LowerBoundary}
+        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,LowerBoundary}}
     end
 
     @testset "boundary" begin

--- a/test/SbpOperators/boundary_conditions/sat_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundary_conditions/sat_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 using Test
 
-using Sbplib.Grids
-using Sbplib.LazyTensors
-using Sbplib.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
+using Diffinitive.SbpOperators
 
 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
 

--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,11 +1,11 @@
 using Test
 
-using Sbplib.LazyTensors
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.RegionIndices
-import Sbplib.SbpOperators.Stencil
-import Sbplib.SbpOperators.BoundaryOperator
+using Diffinitive.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
+import Diffinitive.SbpOperators.Stencil
+import Diffinitive.SbpOperators.BoundaryOperator
 
 
 @testset "BoundaryOperator" begin
@@ -13,12 +13,12 @@
     g_1D = EquidistantGrid(range(0,1,length=11))
 
     @testset "Constructors" begin
-        @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T
-        @test BoundaryOperator(g_1D, closure_stencil, Upper()) isa LazyTensor{T,0,1} where T
+        @test BoundaryOperator(g_1D, closure_stencil, LowerBoundary()) isa LazyTensor{T,0,1} where T
+        @test BoundaryOperator(g_1D, closure_stencil, UpperBoundary()) isa LazyTensor{T,0,1} where T
     end
 
-    op_l = BoundaryOperator(g_1D, closure_stencil, Lower())
-    op_r = BoundaryOperator(g_1D, closure_stencil, Upper())
+    op_l = BoundaryOperator(g_1D, closure_stencil, LowerBoundary())
+    op_r = BoundaryOperator(g_1D, closure_stencil, UpperBoundary())
 
     @testset "Sizes" begin
         @test domain_size(op_l) == (11,)

--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
-using Sbplib.RegionIndices
-using Sbplib.SbpOperators: BoundaryOperator, Stencil
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
+using Diffinitive.RegionIndices
+using Diffinitive.SbpOperators: BoundaryOperator, Stencil
 
 @testset "boundary_restriction" begin
 	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
@@ -14,19 +14,19 @@
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D,stencil_set,Lower())
-            @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Lower())
-            @test e_l isa BoundaryOperator{T,Lower} where T
+            e_l = boundary_restriction(g_1D,stencil_set,LowerBoundary())
+            @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),LowerBoundary())
+            @test e_l isa BoundaryOperator{T,LowerBoundary} where T
             @test e_l isa LazyTensor{T,0,1} where T
 
-            e_r = boundary_restriction(g_1D,stencil_set,Upper())
-            @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Upper())
-            @test e_r isa BoundaryOperator{T,Upper} where T
+            e_r = boundary_restriction(g_1D,stencil_set,UpperBoundary())
+            @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),UpperBoundary())
+            @test e_r isa BoundaryOperator{T,UpperBoundary} where T
             @test e_r isa LazyTensor{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}())
+            e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,UpperBoundary}())
             @test e_w isa InflatedTensor
             @test e_w isa LazyTensor{T,1,2} where T
         end

--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
-using Sbplib.RegionIndices
-import Sbplib.SbpOperators.BoundaryOperator
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
+using Diffinitive.RegionIndices
+import Diffinitive.SbpOperators.BoundaryOperator
 
 @testset "normal_derivative" begin
     g_1D = equidistant_grid(0.0, 1.0, 11)
@@ -12,19 +12,19 @@
     @testset "normal_derivative" begin
     	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            d_l = normal_derivative(g_1D, stencil_set, Lower())
-            @test d_l == normal_derivative(g_1D, stencil_set, Lower())
-            @test d_l isa BoundaryOperator{T,Lower} where T
+            d_l = normal_derivative(g_1D, stencil_set, LowerBoundary())
+            @test d_l == normal_derivative(g_1D, stencil_set, LowerBoundary())
+            @test d_l isa BoundaryOperator{T,LowerBoundary} where T
             @test d_l isa LazyTensor{T,0,1} where T
         end
         @testset "2D" begin
-            d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
-            d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}())
+            d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,UpperBoundary}())
             Ix = IdentityTensor{Float64}((size(g_2D)[1],))
             Iy = IdentityTensor{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower())
-            d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper())
-            @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_2D.grids[1], stencil_set, LowerBoundary())
+            d_r = normal_derivative(g_2D.grids[2], stencil_set, UpperBoundary())
+            @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa LazyTensor{T,1,2} where T

--- a/test/SbpOperators/stencil_set_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/stencil_set_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
 using TOML
-using Sbplib.SbpOperators
+using Diffinitive.SbpOperators
 
-import Sbplib.SbpOperators.Stencil
-import Sbplib.SbpOperators.NestedStencil
+import Diffinitive.SbpOperators.Stencil
+import Diffinitive.SbpOperators.NestedStencil
 
 @testset "readoperator" begin
     toml_str = """

--- a/test/SbpOperators/stencil_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/stencil_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
-using Sbplib.SbpOperators
+using Diffinitive.SbpOperators
 using StaticArrays
-import Sbplib.SbpOperators.Stencil
-import Sbplib.SbpOperators.NestedStencil
-import Sbplib.SbpOperators.scale
-import Sbplib.SbpOperators: apply_stencil, apply_stencil_backwards
+import Diffinitive.SbpOperators.Stencil
+import Diffinitive.SbpOperators.NestedStencil
+import Diffinitive.SbpOperators.scale
+import Diffinitive.SbpOperators: apply_stencil, apply_stencil_backwards
 
 @testset "Stencil" begin
     s = Stencil(-2:2, (1.,2.,2.,3.,4.))

--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,9 +1,9 @@
 using Test
 
-using Sbplib.LazyTensors
-using Sbplib.SbpOperators
-import Sbplib.SbpOperators: ConstantInteriorScalingOperator
-using Sbplib.Grids
+using Diffinitive.LazyTensors
+using Diffinitive.SbpOperators
+import Diffinitive.SbpOperators: ConstantInteriorScalingOperator
+using Diffinitive.Grids
 
 @testset "ConstantInteriorScalingOperator" begin
     @test ConstantInteriorScalingOperator(1, (2,3), 10) isa ConstantInteriorScalingOperator{Int,2}

--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,17 +1,17 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
-using Sbplib.SbpOperators: Stencil
+using Diffinitive.SbpOperators: Stencil
 
-using Sbplib.SbpOperators: dissipation_interior_weights
-using Sbplib.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil
-using Sbplib.SbpOperators: midpoint, midpoint_transpose
-using Sbplib.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size
-using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
-using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
+using Diffinitive.SbpOperators: dissipation_interior_weights
+using Diffinitive.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil
+using Diffinitive.SbpOperators: midpoint, midpoint_transpose
+using Diffinitive.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size
+using Diffinitive.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
+using Diffinitive.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
 
 
 @testset "undivided_skewed04" begin

--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,11 +1,11 @@
 using Test
 
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
-using Sbplib.SbpOperators: closure_size, Stencil, VolumeOperator
+using Diffinitive.SbpOperators: closure_size, Stencil, VolumeOperator
 
 @testset "first_derivative" begin
     @testset "Constructors" begin

--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
-import Sbplib.SbpOperators.VolumeOperator
+import Diffinitive.SbpOperators.VolumeOperator
 
 # TODO: Refactor these test to look more like the tests in first_derivative_test.jl.
 

--- a/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
-using Sbplib.Grids
-using Sbplib.LazyTensors
-using Sbplib.SbpOperators
-using Sbplib.RegionIndices
-using Sbplib.SbpOperators: NestedStencil, CenteredNestedStencil, SecondDerivativeVariable
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.RegionIndices
+using Diffinitive.SbpOperators: NestedStencil, CenteredNestedStencil, SecondDerivativeVariable
 
 using LinearAlgebra
 

--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
-import Sbplib.SbpOperators.ConstantInteriorScalingOperator
+import Diffinitive.SbpOperators.ConstantInteriorScalingOperator
 
 @testset "Diagonal-stencil inner_product" begin
     Lx = π/2.

--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,10 +1,10 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
-import Sbplib.SbpOperators.ConstantInteriorScalingOperator
+import Diffinitive.SbpOperators.ConstantInteriorScalingOperator
 
 @testset "Diagonal-stencil inverse_inner_product" begin
     Lx = π/2.

--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,8 +1,8 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
 @testset "Laplace" begin
     # Default stencils (4th order)

--- a/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,11 +1,11 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
-import Sbplib.SbpOperators.Stencil
-import Sbplib.SbpOperators.StencilOperatorDistinctClosures
+import Diffinitive.SbpOperators.Stencil
+import Diffinitive.SbpOperators.StencilOperatorDistinctClosures
 
 @testset "StencilOperatorDistinctClosures" begin
     g = equidistant_grid(0., 1., 11)

--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,14 +1,14 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.RegionIndices
+using Diffinitive.LazyTensors
 
-import Sbplib.SbpOperators.Stencil
-import Sbplib.SbpOperators.VolumeOperator
-import Sbplib.SbpOperators.odd
-import Sbplib.SbpOperators.even
+import Diffinitive.SbpOperators.Stencil
+import Diffinitive.SbpOperators.VolumeOperator
+import Diffinitive.SbpOperators.odd
+import Diffinitive.SbpOperators.even
 
 
 @testset "VolumeOperator" begin

--- a/test/StaticDicts/StaticDicts_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/StaticDicts/StaticDicts_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -1,5 +1,5 @@
 using Test
-using Sbplib.StaticDicts
+using Diffinitive.StaticDicts
 
 @testset "StaticDicts" begin
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/ext/sparse_array_kit_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,35 @@
+using Test
+
+using Diffinitive
+using Diffinitive.Grids
+using Diffinitive.SbpOperators
+
+using SparseArrayKit
+using Tokens
+using Tullio
+
+
+@testset "SparseArray" begin
+    g = equidistant_grid((0,0),(1,2), 20,30)
+    stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+
+
+    @testset let Δ = laplace(g, stencil_set), M = SparseArray(Δ)
+        @test ndims(M) == 4
+        @test size(M) == (20,30,20,30)
+
+        v = rand(size(g)...)
+        @tullio Mv[i,j] := M[i,j,k,l]*v[k,l]
+
+        @test Mv ≈ Δ*v
+    end
+
+    @testset let dₙ = normal_derivative(g, stencil_set,CartesianBoundary{1,LowerBoundary}()), M = SparseArray(dₙ)
+        @test ndims(M) == 3
+        @test size(M) == (30,20,30)
+
+        v = rand(size(g)...)
+        @tullio Mv[i] := M[i,j,k]*v[j,k]
+        @test Mv ≈ dₙ*v
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/ext/sparse_arrays_test.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -0,0 +1,34 @@
+using Test
+
+using Diffinitive
+using Diffinitive.Grids
+using Diffinitive.SbpOperators
+
+using SparseArrays
+using Tokens
+
+
+@testset "SparseArray" begin
+    g = equidistant_grid((0,0),(1,2), 20,30)
+    stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+
+
+    @testset let Δ = laplace(g, stencil_set), M = sparse(Δ)
+        @test ndims(M) == 2
+        @test size(M) == (20*30,20*30)
+
+        v = rand(size(g)...)
+
+        Mv = M*reshape(v,:)
+        @test Mv ≈ reshape(Δ*v,:)
+    end
+
+    @testset let dₙ = normal_derivative(g, stencil_set,CartesianBoundary{1,LowerBoundary}()), M = sparse(dₙ)
+        @test ndims(M) == 2
+        @test size(M) == (30,20*30)
+
+        v = rand(size(g)...)
+        Mv = M*reshape(v,:)
+        @test Mv ≈ reshape(dₙ*v,:)
+    end
+end

--- a/test/runtests.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/runtests.jl	Wed Sep 11 15:41:58 2024 +0200
@@ -45,7 +45,7 @@
     end
 end
 
-testsetname = isempty(ARGS) ? "Sbplib.jl" : "["*join(ARGS, ", ")*"]"
+testsetname = isempty(ARGS) ? "Diffinitive.jl" : "["*join(ARGS, ", ")*"]"
 
 @testset "$testsetname" begin
     run_testfiles(ARGS)