Mercurial > repos > public > sbplib_julia

changeset 1359:646027afe74b bugfix/lazytensors

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 20 May 2023 14:33:25 +0200
parents 4c0bc52e170f (current diff) 150313ed2cae (diff)
children a6918dfb0cf5
files
diffstat 46 files changed, 1253 insertions(+), 839 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Tue Apr 04 21:46:06 2023 +0200
+++ b/Manifest.toml	Sat May 20 14:33:25 2023 +0200
@@ -2,13 +2,13 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "b024d6898b484706c36ee3b2a041918f3a9d2088"
+project_hash = "a36735c53cfa4453f39635046eeaa47a4ea1231b"
 
 [[deps.Adapt]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c"
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.5.0"
+version = "3.6.1"
 
 [[deps.ArrayInterface]]
 deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"]
@@ -27,9 +27,9 @@
 
 [[deps.Compat]]
 deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04"
+git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.6.0"
+version = "4.6.1"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -77,6 +77,12 @@
 deps = ["SHA", "Serialization"]
 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"
@@ -96,9 +102,24 @@
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c"
+git-tree-sha1 = "08be5ee09a7632c32695d954a602df96a877bf0d"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.3"
+version = "0.8.6"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
+git-tree-sha1 = "7756ce473bd10b67245bdebdc8d8670a85f6230b"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.5.18"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.0"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
--- a/Notes.md	Tue Apr 04 21:46:06 2023 +0200
+++ b/Notes.md	Sat May 20 14:33:25 2023 +0200
@@ -141,13 +141,13 @@
  - [ ] 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.
- - [ ] How do we handle mixes of periodic and non-periodic grids? Seems it should be supported on the grid level and on the 1d operator level. Between there it should be transparent.
  - [ ] 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 probabily be included in the grid module. This allows new grid types to come with their own 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.
@@ -198,80 +198,83 @@
 
 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.
 
-## Vector valued grid functions
-Från slack konversation:
+## 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.
+
+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.
 
-Jonatan Werpers:
-Med vektorvärda gridfunktioner vill vi ju fortfarande att grid funktionen ska vara till exempel AbstractArray{LitenVektor,2}
-Och att man ska kunna göra allt man vill med LitenVektor
-typ addera, jämföra osv
-Och då borde points returnera AbstractArray{LitenVektor{Float,2},2} för ett 2d nät
-Men det kanske bara ska vara Static arrays?
+* 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.
 
-Vidar Stiernström:
-Ja, jag vet inte riktigt vad som är en rimlig representation
-Du menar en vektor av static arrays då?
+Below is a partial implementation of `lmap` with some ideas
+```julia
+struct LazyMapping{T,IT,F}
+    f::F
+    indexable_iterator::IT # ___
+end
 
-Jonatan Werpers:
-Ja, att LitenVektor är en StaticArray
+function LazyMapping(f,I)
+    IT = eltype(I)
+    T = f(zero(T))
+    F = typeof(f)
 
-Vidar Stiernström:
-Tuplar känns typ rätt inuitivt för att representera värdet i en punkt
-men
-det suger att man inte har + och - för dem
+    return LazyMapping{T,IT,F}(f,I)
+end
 
-Jonatan Werpers:
-Ja precis
+getindex(lm::LazyMapping, I...) = lm.f(lm.I[I...])
+# indexabl interface
+# iterable has shape
 
-Vidar Stiernström:
-så kanske är bra med static arrays i detta fall
+iterate(lm::LazyMapping) = _lazy_mapping_iterate(lm, iterate(lm.I))
+iterate(lm::LazyMapping, state) = _lazy_mapping_iterate(lm, iterate(lm.I, state))
 
-Jonatan Werpers:
-Man vill ju kunna köra en Operator rakt på och vara klar eller?
+_lazy_mapping_iterate(lm, ::Nothing) = nothing
+_lazy_mapping_iterate(lm, (next, state)) = lm.f(next), state
 
-Vidar Stiernström:
-Har inte alls tänkt på hur det vi gör funkar mot vektorvärda funktioner
-men känns som staticarrays är hur man vill göra det
-tuplar är ju immutable också
-blir jobbigt om man bara agerar på en komponent då
+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.
+
+## Multiblock implementation
+We want multiblock things to work very similarly to regular one block things.
 
-Jonatan Werpers:
-Hm…
-Tål att tänkas på
-Men det lär ju bli mer indirektion med mutables eller?
-Hur fungerar det?
-Det finns ju hur som helst både SVector och MVector i StaticArrays
+### Grid functions
+Should probably support a nested indexing so that we first have an index for subgrid and then an index for nodes on that grid. E.g `g[1,2][2,3]` or `g[3][43,21]`.
+
+We could also possibly provide a combined indexing style `g[1,2,3,4]` where the first group of indices are for the subgrid and the remaining are for the nodes.
 
-Vidar Stiernström:
-När vi jobbat i c/c++ och kollat runt lite hur man brukar göra så lagrar man i princip alla sina obekanta i en lång vektor och så får man specificera i funktioerna vilken komponent man agerar på och till vilken man skriver
-så man lagrar grejer enl: w = [u1, v1, u2, v2, …] i 1D.
-Men alltså har ingen aning hur julia hanterar detta
+We should make sure the underlying buffer for gridfunctions are continuously stored and are easy to convert to, so that interaction with for example DifferentialEquations is simple and without much boilerplate.
+
+#### `map` and `collect` and nested indexing
+We need to make sure `collect`, `map` and a potential lazy map work correctly through the nested indexing.
 
-Jonatan Werpers:
-Det vi är ute efter kanske är att en grid funcktion är en AbstractArray{T,2} where T<:NågotSomViKanRäknaMed
-Och så får den typen var lite vad som helst.
+### Tensor applications
+Should behave as grid functions
 
-Vidar Stiernström:
-Tror det kan vara farligt att ha nåt som är AbstractArray{LitenArray{NDof},Dim}
-Jag gissar att det kompilatorn vill ha är en stor array med doubles
+### LazyTensors
+Could be built as a tuple or array of LazyTensors for each grid with a simple apply function.
+
+Nested indexing for these is problably not needed unless it simplifies their own implementation.
 
-Jonatan Werpers:
-Och sen är det upp till den som använder grejerna att vara smart
-Vill man vara trixig kan man väl då imlementera SuperHaxxorGridFunction <: AbstractArray{Array{…},2} som lagrar allt linjärt eller något sånt
-Det kommer väl lösa sig när man börjar implementera vektorvärda saker
-Euler nästa!
-New
-Vidar Stiernström:
-Det vore skönt att inte behöva skriva såhär varje gång man testar mot en tupel :smile: @test [gp[i]...] ≈ [p[i]...] atol=5e-13
+Possibly useful to provide a simple type that doesn't know about connections between the grids. Antother type can include knowledge of the.
 
-Jonatan Werpers:
-https://github.com/JuliaArrays/ArraysOfArrays.jl
-https://github.com/jw3126/Setfield.jl
+We have at least two option for how to implement them:
+* Matrix of LazyTensors
+* Looking at the grid and determining what the apply should do.
+
+### Overall design implications of nested indices
+If some grids accept nested indexing there might be a clash with how LazyArrays work. It would be nice if the grid functions and lazy arrays that actually are arrays can be AbstractArray and things can be relaxed for nested index types.
+
+## Vector valued grid functions
 
 ### Test-applikationer
-div och grad operationer
+div- och grad-operationer
 
-Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensor fält och `div()` på sista.
+Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensorfält och `div()` på sista.
 Om man generaliserar kanske `∇` i så fall bara lägger till en dimension i början.
 
 Kan vi implementera `⋅`(\cdot) så att de fungerar som man vill för både tensor-fält och tensor-operatorer?
@@ -279,8 +282,8 @@
 Är `∇` ett tensor-fält av tensor-operatorer? Vad är ett tensor-fält i vår kod? Är det en special-fall av en tensor-mapping?
 
 ### Grid-funktionen
-Grid-funktionon 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.
+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.
 
@@ -292,8 +295,6 @@
 gf[2,3][2] # x̄[2] för en viss gridpunkt
 ```
 
-Note: Behöver bestämma om eval on skickar in `x̄` eller `x̄...` till `f`. Eller om man kan stödja båda.
-
 ### 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.
@@ -309,63 +310,35 @@
 
 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.
 
-### Ratade alternativ:
-
-
-#### 2.AbstractArray{T,2+1} where T (NOPE!)
-Blir inte den här. Bryter mot alla tankar om hur grid funktioner ska fungera. Om de tillåts ha en annan dimension än nätet blir allt hemskt.
-
-Man låter helt enkelt arrayen ha en extra dimension. En fördel är att man har en väldigt "native" typ. En nackdel kan vara att det eventuellt blir rörigt vilken dimension olika operatorer ska agera på. I värsta fall behöver vi "kroneckra in" de tillagda dimensionerna. Vektorfältets index kommer också att bli det första eftersom vi vill att de ska lagras kontinuerligt i minnet pga chachen. (Går kanske att lösa med en custom typ men då krånglar man till det för sig). En fördel skulle vara att man enkelt får ut olika komponenter.
-
-Syntax:
-```
-gf = eval_on_grid(g,f)
-gf[:,2,3] # Hela vektorn för en gridpunkt
-gf[2,2,3] # Andra komponenten av vektor fältet i en punkt.
-gf[2,:,:] #
-```
-
-### Evaluering av funktioner på nät
-Hur ska man skriva funktioner som evalueras på nätet? `f(x,y) = ...` eller `f(x̄) = ...`? Eller båda? Kan eval_on_grid se skillnad eller får användaren specificera?
-
-```
-f(x,y) = [x^2, y^2]
-f(x̄) = [x̄[1]^2, x̄[2]^2]
-```
-
-Påverkas detta av hur vi förväntar oss kunna skapa lata gridfunktioner?
-
 ### Komponenter som gridfunktioner
-En viktig operation för vektor fält är att kunna få ut komponenter som grid-funktioner. Detta behöver antagligen kunna ske lazy.
+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.
 
-Skulle vara en fördel om det är hyffsat generiskt så att en eventuell användare kan utöka det enkelt om de har någon egen exotisk typ. Eller ska man vila helt på
 
-Syntax:
-```
-gf = eval(...)
-component(gf,2) # Andra komponenten av en vektor
-component(gf,2,3) # (2,3) elementet av en matris
-component(gf,:,2) # Andra kolumnen av en matris
-@ourview gf[:,:][2]
-```
+### 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
+u = [ [u1[x1], u2[x1]] , [u1[x2], u2[x2]], ... [u1[xN], u2[xN]]]. Detta brukar kallas array of structs (AoS). Alternativet är struct of arrays (SoA), där man har alla gridpunkter för en given komponent u = [[u1[x1], u1[x2]],... u1[xN]], [u2[x1], u2[x2], ... u2[xN]]].
+
+Personligen tycker jag att AoS känns som den mer naturliga representationen? Det skulle göra det enklarare att parallelisera en vektorvärd gridfunktion över gridpunkterna, och om man opererar på olika komponenter i samma funktion så är det också bra ur en minnesaccess-synpunkt då dessa kommer ligga nära vandra i minnet. Problemet är att AoS sabbar vektorisering på CPU då två gridpunkter i en komponent ligger långt bort från varandra. Efter lite eftersökningar (och efter att snackat lite med Ossian) så verkar det ändå som att AoS är dåligt på GPU, där man vill att trådar typiskt sett utföra samma operation på närliggande minne.
+
+Vad tänker du kring detta ur ett interface-perspektiv? Jag hittade paketet  https://github.com/JuliaArrays/StructArrays.jl som verkar erbjuda AoS-interface men SoA-minneslayout så det kanske kan vara något vi kan använda? Inte native-stödd på samma sätt som SVector, men verkar iaf utvecklas aktivt.
 
-## Grids embedded in higher dimensions
+[Efter telefonsamtal] För optimal prestanda behöver vi antagligen se till att man kan räkna ut varje komponent i en punkt individuellt. Detta så att man har frihet att till exempel låta den innersta loopen hålla komponentindexet konstant för att underlätta intruktionsvektorisering.
+
 
-For grids generated by asking for boundary grids for a regular grid, it would
-make sense if these grids knew they were embedded in a higher dimension. They
-would return coordinates in the full room. This would make sense when
-drawing points for example, or when evaluating functions on the boundary.
+[Vidare tankar]
+ * Det borde bara vara output-gridfunktionen som behöver special-indexeras? Det viktiga på inputsidan är att den är lagrad på rätt sätt i minnet.
+ * Det borde inte vara några problem att behålla det "optimala" interfacet (gf[1,1,1][2]) till gridfunktionerna. Om man verkligen behöver kan skapa parallella indexeringsmetoder som gör det man behöver, i.e, "deep indexing".
+ * Det är inte säkert att vi behöver göra något speciellt på outputsidan överhuvudtaget. Det känns inte orimligt att kompilatorn skulle kunna optimera bort den koden som räknar ut onödiga komponenter.
+ * Om vi behöver special-indexering kommer till exempel LazyTensorApplication att behöva implementera det.
+ * För att komma vidare med något mer avancerat behöver vi antagligen implementera några operatorer som ger och agerar på vektorvärda funktioner. Tex grad, elastiska operatorn, andra?
 
-Implementation of this is an issue that requires some thought. Adding an extra
-"Embedded" type for each grid would make it easy to understand each type but
-contribute to "type bloat". On the other hand adapting existing types to
-handle embeddedness would complicate the now very simple grid types. Are there
-other ways of doing the implentation?
 
 ## Performance measuring
 We should be measuring performance early. How does our effective cpu and memory bandwidth utilization compare to peak performance?
--- a/Project.toml	Tue Apr 04 21:46:06 2023 +0200
+++ b/Project.toml	Sat May 20 14:33:25 2023 +0200
@@ -4,6 +4,7 @@
 version = "0.1.0"
 
 [deps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
--- a/TODO.md	Tue Apr 04 21:46:06 2023 +0200
+++ b/TODO.md	Sat May 20 14:33:25 2023 +0200
@@ -1,5 +1,7 @@
 # TODO
 
+## Organization
+ - [ ] Split up Notes.md in several files
 
 ## Coding
  - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
@@ -26,3 +28,16 @@
  - [ ] Kolla att vi gör boundschecks överallt och att de är markerade med @boundscheck
  - [ ] Kolla att vi har @inline på rätt ställen
  - [ ] Profilera
+
+
+### Grids
+
+ - [ ] Multiblock grids
+ - [ ] Periodic grids
+ - [ ] Grids with modified boundary closures
+
+
+### Benchmarks
+ - [ ] Benchmarks for all grid indexing (focused on allocation)
+ - [ ] Benchmarks for indexing of lazy grid functions
+ - [ ] Add benchmarks for range type in EquidistantGrid. (LinRange vs StepRange)
--- a/docs/Manifest.toml	Tue Apr 04 21:46:06 2023 +0200
+++ b/docs/Manifest.toml	Sat May 20 14:33:25 2023 +0200
@@ -10,10 +10,10 @@
 version = "0.0.1"
 
 [[deps.Adapt]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c"
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.5.0"
+version = "3.6.1"
 
 [[deps.ArrayInterface]]
 deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"]
@@ -35,9 +35,9 @@
 
 [[deps.Compat]]
 deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04"
+git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.6.0"
+version = "4.6.1"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -119,9 +119,9 @@
 
 [[deps.Parsers]]
 deps = ["Dates", "SnoopPrecompile"]
-git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b"
+git-tree-sha1 = "478ac6c952fddd4399e71d4779797c538d0ff2bf"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.5.4"
+version = "2.5.8"
 
 [[deps.Preferences]]
 deps = ["TOML"]
@@ -141,12 +141,18 @@
 deps = ["SHA", "Serialization"]
 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.Sbplib]]
-deps = ["TOML", "TiledIteration"]
+deps = ["StaticArrays", "TOML", "TiledIteration"]
 path = ".."
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
 version = "0.1.0"
@@ -169,9 +175,24 @@
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c"
+git-tree-sha1 = "08be5ee09a7632c32695d954a602df96a877bf0d"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.3"
+version = "0.8.6"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
+git-tree-sha1 = "7756ce473bd10b67245bdebdc8d8670a85f6230b"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.5.18"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.0"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
--- a/docs/make.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/docs/make.jl	Sat May 20 14:33:25 2023 +0200
@@ -26,6 +26,7 @@
 pages = [
     "Home" => "index.md",
     "operator_file_format.md",
+    "grids_and_grid_functions.md",
     "Submodules" => [
         "submodules/grids.md",
         "submodules/diff_ops.md",
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/grids_and_grid_functions.md	Sat May 20 14:33:25 2023 +0200
@@ -0,0 +1,17 @@
+# Grids and grid functions
+
+The submodule `Grids` aims to provide types and logic for all types of grids that are useful for implementing summation-by-parts difference methods. It provides an abstract top level type `Grid` which defines a broad interface for how a general grid is supposed to work. Currently only equidistant grids are supported, but the basic structure supports implementations of curvilinear grids, multi-block grids, periodic grids and much more.
+
+The module also has functionality for creating and working with grid functions.
+
+## Interface for grids
+All grids are expected to work as a grid function for the coordinate function, and thus implements Julia's Indexing- and Iteration-interfaces. Notably they are *not* abstract arrays because that inteface is too restrictive for the types of grids we wish to implement.
+
+## To write about
+<!-- # TODO: -->
+* Grid functions
+  * Basic structure
+     * Indexing
+  * Curvilinear
+  * Multiblock
+  * Vector valued grid functions
--- a/src/Grids/Grids.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/Grids/Grids.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,31 +1,45 @@
 module Grids
 
 using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+using StaticArrays
 
 # Grid
 export Grid
-export dims
-export points
-export evalOn
+export coordinate_size
+export component_type
+
+export TensorGrid
+export ZeroDimGrid
+
+export TensorGridBoundary
+
+export grid_id
+export boundary_id
+
+export eval_on
+export getcomponent
 
 # BoundaryIdentifier
 export BoundaryIdentifier
-export CartesianBoundary
-export dim
-export region
+
 
 # EquidistantGrid
 export EquidistantGrid
 export spacing
 export inverse_spacing
-export restrict
 export boundary_identifiers
 export boundary_grid
 export refine
 export coarsen
+export equidistant_grid
+export CartesianBoundary
+
+abstract type BoundaryIdentifier end
 
 include("grid.jl")
-include("boundary_identifier.jl")
+include("tensor_grid.jl")
 include("equidistant_grid.jl")
+include("zero_dim_grid.jl")
 
 end # module
--- a/src/Grids/boundary_identifier.jl	Tue Apr 04 21:46:06 2023 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-
-abstract type BoundaryIdentifier end
-
-struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end
-dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim
-region(::CartesianBoundary{Dim, R}) where {Dim, R} = R()
\ No newline at end of file
--- a/src/Grids/equidistant_grid.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/Grids/equidistant_grid.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,71 +1,32 @@
-
 """
-    EquidistantGrid{Dim,T<:Real} <: Grid
-
-`Dim`-dimensional equidistant grid with coordinates of type `T`.
-"""
-struct EquidistantGrid{Dim,T<:Real} <: Grid
-    size::NTuple{Dim, Int}
-    limit_lower::NTuple{Dim, T}
-    limit_upper::NTuple{Dim, T}
+    EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1}
 
-    function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T}
-        if any(size .<= 0)
-            throw(DomainError("all components of size must be postive"))
-        end
-        if any(limit_upper.-limit_lower .<= 0)
-            throw(DomainError("all side lengths must be postive"))
-        end
-        return new{Dim,T}(size, limit_lower, limit_upper)
-    end
-end
+A one-dimensional equidistant grid. Most users are expected to use
+[`equidistant_grid`](@ref) for constructing equidistant grids.
+
+See also: [`equidistant_grid`](@ref)
 
 
+## Note
+The type of range used for the points can likely impact performance.
 """
-    EquidistantGrid(size, limit_lower, limit_upper)
-
-Construct an equidistant grid with corners at the coordinates `limit_lower` and
-`limit_upper`.
-
-The length of the domain sides are given by the components of
-`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and `limit_upper=(1,2)` the domain is defined
-as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative.
-
-The number of equidistantly spaced points in each coordinate direction are given
-by the tuple `size`.
-"""
-function EquidistantGrid(size, limit_lower, limit_upper)
-    return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper)
+struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1}
+    points::R
 end
 
-
-"""
-    EquidistantGrid{T}()
-
-Constructs a 0-dimensional grid.
-"""
-EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid
-
-
-"""
-    EquidistantGrid(size::Int, limit_lower::T, limit_upper::T)
+# Indexing interface
+Base.getindex(g::EquidistantGrid, i) = g.points[i]
+Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
+Base.firstindex(g::EquidistantGrid) = firstindex(g.points)
+Base.lastindex(g::EquidistantGrid) = lastindex(g.points)
 
-Convenience constructor for 1D grids.
-"""
-function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T
-	return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
-end
-
-Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
+# Iteration interface
+Base.iterate(g::EquidistantGrid) = iterate(g.points)
+Base.iterate(g::EquidistantGrid, state) = iterate(g.points, state)
 
-Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
-
-Base.size(g::EquidistantGrid) = g.size
-
-Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
-
-
-
+Base.IteratorSize(::Type{<:EquidistantGrid}) = Base.HasShape{1}()
+Base.length(g::EquidistantGrid) = length(g.points)
+Base.size(g::EquidistantGrid) = size(g.points)
 
 
 """
@@ -73,7 +34,7 @@
 
 The spacing between grid points.
 """
-spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1)
+spacing(g::EquidistantGrid) = step(g.points)
 
 
 """
@@ -81,111 +42,98 @@
 
 The reciprocal of the spacing between grid points.
 """
-inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid)
+inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
 
 
-"""
-    points(grid::EquidistantGrid)
-
-The point of the grid as an array of tuples with the same dimension as the grid.
-The points are stored as [(x1,y1), (x1,y2), … (x1,yn);
-						  (x2,y1), (x2,y2), … (x2,yn);
-						  	⋮		 ⋮            ⋮
-						  (xm,y1), (xm,y2), … (xm,yn)]
-"""
-function points(grid::EquidistantGrid)
-    indices = Tuple.(CartesianIndices(grid.size))
-    h = spacing(grid)
-    return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
-end
+boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
+boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
+boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
 
 
 """
-    restrict(::EquidistantGrid, dim)
+    refine(g::EquidistantGrid, r::Int)
 
-Pick out given dimensions from the grid and return a grid for them.
+The grid where `g` is refined by the factor `r`. The factor is applied to the number of
+intervals, i.e., 1 less than the size of  `g`.
+
+See also: [`coarsen`](@ref)
 """
-function restrict(grid::EquidistantGrid, dim)
-    size = grid.size[dim]
-    limit_lower = grid.limit_lower[dim]
-    limit_upper = grid.limit_upper[dim]
-
-    return EquidistantGrid(size, limit_lower, limit_upper)
+function refine(g::EquidistantGrid, r::Int)
+    new_sz = (length(g) - 1)*r + 1
+    return EquidistantGrid(change_length(g.points, new_sz))
 end
 
+"""
+    coarsen(g::EquidistantGrid, r::Int)
 
-"""
-    orthogonal_dims(grid::EquidistantGrid,dim)
+The grid where `g` is coarsened by the factor `r`. The factor is applied to the number of
+intervals, i.e., 1 less than the size of `g`. If the number of
+intervals are not divisible by `r` an error is raised.
 
-Returns the dimensions of grid orthogonal to that of dim.
+See also: [`refine`](@ref)
 """
-function orthogonal_dims(grid::EquidistantGrid, dim)
-    orth_dims = filter(i -> i != dim, dims(grid))
-	if orth_dims == dims(grid)
-		throw(DomainError(string("dimension ",string(dim)," not matching grid")))
-	end
-    return orth_dims
+function coarsen(g::EquidistantGrid, r::Int)
+    if (length(g)-1)%r != 0
+        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
+    end
+
+    new_sz = (length(g) - 1)÷r + 1
+
+    return EquidistantGrid(change_length(g.points, new_sz))
 end
 
 
 """
-    boundary_identifiers(::EquidistantGrid)
+    equidistant_grid(size::Dims, limit_lower, limit_upper)
+
+Construct an equidistant grid with corners at the coordinates `limit_lower` and
+`limit_upper`.
+
+The length of the domain sides are given by the components of
+`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and
+`limit_upper=(1,2)` the domain is defined as `(-1,1)x(0,2)`. The side lengths
+of the grid are not allowed to be negative.
 
-Returns a tuple containing the boundary identifiers for the grid, stored as
-	(CartesianBoundary(1,Lower),
-	 CartesianBoundary(1,Upper),
-	 CartesianBoundary(2,Lower),
-	 ...)
+The number of equispaced points in each coordinate direction are given
+by the tuple `size`.
+
+Note: If `limit_lower` and `limit_upper` are integers and `size` would allow a
+completely integer grid, `equidistant_grid` will still return a floating point
+grid. This simlifies the implementation and avoids certain surprise
+behaviours.
 """
-boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
-
+function equidistant_grid(size::Dims, limit_lower, limit_upper)
+    gs = map(equidistant_grid, size, limit_lower, limit_upper)
+    return TensorGrid(gs...)
+end
 
 """
-    boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
+    equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
 
-Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
-orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
-grid is a zero-dimensional grid.
+Constructs a 1D equidistant grid.
 """
-function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
-    orth_dims = orthogonal_dims(grid, dim(id))
-    return restrict(grid, orth_dims)
+function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
+    if any(size .<= 0)
+        throw(DomainError("size must be postive"))
+    end
+
+    if any(limit_upper.-limit_lower .<= 0)
+        throw(DomainError("side length must be postive"))
+    end
+	return EquidistantGrid(range(limit_lower, limit_upper, length=size)) # TBD: Should it use LinRange instead?
 end
-boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
+
+CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
 
 
 """
-    refine(grid::EquidistantGrid, r::Int)
-
-Refines `grid` by a factor `r`. The factor is applied to the number of
-intervals which is 1 less than the size of the grid.
-
-See also: [`coarsen`](@ref)
-"""
-function refine(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
-    new_sz = (sz .- 1).*r .+ 1
-    return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end
-
-
-"""
-    coarsen(grid::EquidistantGrid, r::Int)
+    change_length(r::AbstractRange, n)
 
-Coarsens `grid` by a factor `r`. The factor is applied to the number of
-intervals which is 1 less than the size of the grid. If the number of
-intervals are not divisible by `r` an error is raised.
-
-See also: [`refine`](@ref)
+Change the length of `r` to `n`, keeping the same start and stop.
 """
-function coarsen(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
+function change_length end
 
-    if !all(n -> (n % r == 0), sz.-1)
-        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
-    end
-
-    new_sz = (sz .- 1).÷r .+ 1
-
-    return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end
+change_length(r::UnitRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
+change_length(r::StepRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
+change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
+change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)
--- a/src/Grids/grid.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/Grids/grid.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,27 +1,93 @@
 """
-     Grid
+     Grid{T,D}
+
+A grid with coordinates of type `T`, e.g. `SVector{3,Float64}`, and dimension
+`D`. The grid can be embedded in a higher dimension in which case the number
+of indices and the number of components of the coordinate vectors will be
+different.
+
+All grids are expected to behave as a grid function for the coordinates.
 
-Should implement
-    Base.ndims(grid::Grid)
-    points(grid::Grid)
+`Grids` is top level abstract type for grids. A grid should implement Julia's interfaces for
+indexing and iteration.
+
+## Note
+
+Importantly a grid does not have to be an `AbstractArray`. The reason is to
+allow flexible handling of special types of grids like multi-block grids, or
+grids with special indexing.
+"""
+abstract type Grid{T,D} end
+
+Base.ndims(::Grid{T,D}) where {T,D} = D
+Base.eltype(::Type{<:Grid{T}}) where T = T
 
 """
-abstract type Grid end
-function points end
+    coordinate_size(g)
+
+The lenght of the coordinate vector of `Grid` `g`.
+"""
+coordinate_size(::Type{<:Grid{T}}) where T = _ncomponents(T)
+coordinate_size(g::Grid) = coordinate_size(typeof(g)) # TBD: Name of this function?!
+
+"""
+    component_type(g)
+
+The type of the components of the coordinate vector of `Grid` `g`.
+"""
+component_type(::Type{<:Grid{T}}) where T = eltype(T)
+component_type(g::Grid) = component_type(typeof(g))
 
 """
-    dims(grid::Grid)
+    refine(g::Grid, r)
+
+The grid where `g` is refined by the factor `r`.
+
+See also: [`coarsen`](@ref).
+"""
+function refine end
+
+"""
+    coarsen(g::Grid, r)
 
-A range containing the dimensions of `grid`
+The grid where `g` is coarsened by the factor `r`.
+
+See also: [`refine`](@ref).
 """
-dims(grid::Grid) = 1:ndims(grid)
+function coarsen end
+
+"""
+    boundary_identifiers(g::Grid)
+
+Identifiers for all the boundaries of `g`.
+"""
+function boundary_identifiers end
 
 """
-    evalOn(grid::Grid, f::Function)
+    boundary_grid(g::Grid, id::BoundaryIdentifier)
 
-Evaluate function `f` on `grid`
+The grid for the boundary specified by `id`.
+"""
+function boundary_grid end
+# TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff?
+
 """
-function evalOn(grid::Grid, f::Function)
-    F(x) = f(x...)
-    return F.(points(grid))
+    eval_on(g::Grid, f)
+
+Lazy evaluation `f` on the grid. `f` can either be on the form `f(x,y,...)`
+with each coordinate as an argument, or on the form `f(x̄)` taking a
+coordinate vector.
+
+For concrete array grid functions `map(f,g)` can be used instead.
+"""
+eval_on(g::Grid, f) = eval_on(g, f, Base.IteratorSize(g))
+function eval_on(g::Grid, f, ::Base.HasShape)
+    if hasmethod(f, (Any,))
+        return LazyTensors.LazyFunctionArray((I...)->f(g[I...]), size(g))
+    else
+        return LazyTensors.LazyFunctionArray((I...)->f(g[I...]...), size(g))
+    end
 end
+
+_ncomponents(::Type{<:Number}) = 1
+_ncomponents(T::Type{<:SVector}) = length(T)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/tensor_grid.jl	Sat May 20 14:33:25 2023 +0200
@@ -0,0 +1,99 @@
+"""
+    TensorGrid{T,D} <: Grid{T,D}
+
+A grid constructed as the tensor product of other grids.
+
+Currently only supports grids with the `HasShape`-trait.
+"""
+struct TensorGrid{T,D,GT<:NTuple{N,Grid} where N} <: Grid{T,D}
+    grids::GT
+
+    function TensorGrid(gs...)
+        T = mapreduce(eltype, combined_coordinate_vector_type, gs)
+        D = sum(ndims, gs)
+
+        return new{T,D,typeof(gs)}(gs)
+    end
+end
+
+# Indexing interface
+function Base.getindex(g::TensorGrid, I...)
+    szs = ndims.(g.grids)
+
+    Is = LazyTensors.split_tuple(I, szs)
+    ps = map((g,I)->SVector(g[I...]), g.grids, Is)
+
+    return vcat(ps...)
+end
+
+Base.getindex(g::TensorGrid, I::CartesianIndex) = g[Tuple(I)...]
+
+function Base.eachindex(g::TensorGrid)
+    szs = LazyTensors.concatenate_tuples(size.(g.grids)...)
+    return CartesianIndices(szs)
+end
+
+# Iteration interface
+Base.iterate(g::TensorGrid) = iterate(Iterators.product(g.grids...)) |> _iterate_combine_coords
+Base.iterate(g::TensorGrid, state) = iterate(Iterators.product(g.grids...), state) |> _iterate_combine_coords
+_iterate_combine_coords(::Nothing) = nothing
+_iterate_combine_coords((next,state)) = combine_coordinates(next...), state
+
+Base.IteratorSize(::Type{<:TensorGrid{<:Any, D}}) where D = Base.HasShape{D}()
+Base.eltype(::Type{<:TensorGrid{T}}) where T = T
+Base.length(g::TensorGrid) = sum(length, g.grids)
+Base.size(g::TensorGrid) = LazyTensors.concatenate_tuples(size.(g.grids)...)
+
+
+refine(g::TensorGrid, r::Int) = mapreduce(g->refine(g,r), TensorGrid, g.grids)
+coarsen(g::TensorGrid, r::Int) = mapreduce(g->coarsen(g,r), TensorGrid, g.grids)
+
+"""
+    TensorGridBoundary{N, BID} <: BoundaryIdentifier
+
+A boundary identifier for a tensor grid. `N` Specifies which grid in the
+tensor product and `BID` which boundary on that grid.
+"""
+struct TensorGridBoundary{N, BID} <: BoundaryIdentifier end
+grid_id(::TensorGridBoundary{N, BID}) where {N, BID} = N
+boundary_id(::TensorGridBoundary{N, BID}) where {N, BID} = BID()
+
+"""
+    boundary_identifiers(g::TensorGrid)
+
+Returns a tuple containing the boundary identifiers of `g`.
+"""
+function boundary_identifiers(g::TensorGrid)
+    per_grid = map(eachindex(g.grids)) do i
+        return map(bid -> TensorGridBoundary{i, typeof(bid)}(), boundary_identifiers(g.grids[i]))
+    end
+    return LazyTensors.concatenate_tuples(per_grid...)
+end
+
+
+"""
+    boundary_grid(g::TensorGrid, id::TensorGridBoundary)
+
+The grid for the boundary of `g` specified by `id`.
+"""
+function boundary_grid(g::TensorGrid, id::TensorGridBoundary)
+    local_boundary_grid = boundary_grid(g.grids[grid_id(id)], boundary_id(id))
+    new_grids = Base.setindex(g.grids, local_boundary_grid, grid_id(id))
+    return TensorGrid(new_grids...)
+end
+
+
+function combined_coordinate_vector_type(coordinate_types...)
+    combined_coord_length = mapreduce(_ncomponents, +, coordinate_types)
+    combined_coord_type = mapreduce(eltype, promote_type, coordinate_types)
+
+    if combined_coord_length == 1
+        return combined_coord_type
+    else
+        return SVector{combined_coord_length, combined_coord_type}
+    end
+end
+
+function combine_coordinates(coords...)
+    return mapreduce(SVector, vcat, coords)
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/zero_dim_grid.jl	Sat May 20 14:33:25 2023 +0200
@@ -0,0 +1,27 @@
+"""
+    ZeroDimGrid{T} <: Grid{T,0}
+
+A zero dimensional grid consisting of a single point.
+"""
+struct ZeroDimGrid{T} <: Grid{T,0}
+    point::T
+end
+
+# Indexing interface
+Base.getindex(g::ZeroDimGrid) = g.point
+Base.eachindex(g::ZeroDimGrid) = CartesianIndices(())
+
+# Iteration interface
+Base.iterate(g::ZeroDimGrid) = (g.point, nothing)
+Base.iterate(g::ZeroDimGrid, ::Any) = nothing
+
+Base.IteratorSize(::Type{<:ZeroDimGrid}) = Base.HasShape{0}()
+Base.length(g::ZeroDimGrid) = 1
+Base.size(g::ZeroDimGrid) = ()
+
+
+refine(g::ZeroDimGrid, ::Int) = g
+coarsen(g::ZeroDimGrid, ::Int) = g
+
+boundary_identifiers(g::ZeroDimGrid) = ()
+boundary_grid(g::ZeroDimGrid, ::Any) = throw(ArgumentError("ZeroDimGrid has no boundaries"))
--- a/src/LazyTensors/LazyTensors.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/LazyTensors/LazyTensors.jl	Sat May 20 14:33:25 2023 +0200
@@ -30,6 +30,7 @@
 
 # Composing lazy tensors
 Base.:∘(s::LazyTensor, t::LazyTensor) = TensorComposition(s,t)
+Base.:∘(s::TensorComposition, t::LazyTensor) = s.t1∘(s.t2∘t)
 
 # Outer products of tensors
 ⊗(a::LazyTensor, b::LazyTensor) = LazyOuterProduct(a,b)
--- a/src/LazyTensors/lazy_array.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/LazyTensors/lazy_array.jl	Sat May 20 14:33:25 2023 +0200
@@ -36,7 +36,7 @@
 
 function Base.getindex(lfa::LazyFunctionArray{F,T,D}, I::Vararg{Int,D}) where {F,T,D}
     @boundscheck checkbounds(lfa, I...)
-    return @inbounds lfa.f(I...)
+    return @inbounds @inline lfa.f(I...)
 end
 
 
--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Sat May 20 14:33:25 2023 +0200
@@ -3,9 +3,10 @@
 
 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 prolongation of a zero dimensional array to the whole grid using the same closure stencil.
+`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
+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}
     stencil::Stencil{T,N}
@@ -13,12 +14,13 @@
 end
 
 """
-    BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region)
+    BoundaryOperator(grid::EquidistantGrid, closure_stencil, region)
 
-Constructs the BoundaryOperator with stencil `closure_stencil` for a one-dimensional `grid`, restricting to
-to the boundary specified by `region`.
+Constructs the BoundaryOperator with stencil `closure_stencil` for a
+`EquidistantGrid` `grid`, restricting to to the boundary specified by
+`region`.
 """
-function BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil::Stencil{T,N}, region::Region) where {T,N}
+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])
 end
 
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,25 +1,27 @@
 """
-    boundary_restriction(grid, closure_stencil::Stencil, boundary)
+    boundary_restriction(g, stencil_set::StencilSet, boundary)
+    boundary_restriction(g::TensorGrid, stencil_set::StencilSet, boundary::TensorGridBoundary)
+    boundary_restriction(g::EquidistantGrid, stencil_set::StencilSet, boundary)
 
 Creates boundary restriction operators `e` as `LazyTensor`s on `boundary`
 
-`e` is the restriction of a grid function to `boundary` using a `Stencil` `closure_stencil`.
-`e'` is the prolongation of a grid function on `boundary` to the whole grid using the same `closure_stencil`.
-On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
-a `BoundaryOperator`.
+`e` restricts a grid function on `g` to `boundary` using the 'e' stencil
+in `stencil_set`. `e'` prolongates a grid function on
+`boundary` to the whole grid using the same stencil. On a one-dimensional
+grid, `e` is a `BoundaryOperator`. On a multi-dimensional grid, `e` is the
+inflation of a `BoundaryOperator`.
 
 See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
-function boundary_restriction(grid, closure_stencil, boundary)
-    converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil)
+function boundary_restriction end
 
-    op = BoundaryOperator(restrict(grid, dim(boundary)), converted_stencil, region(boundary))
-    return LazyTensors.inflate(op, size(grid), dim(boundary))
+function boundary_restriction(g::TensorGrid, stencil_set::StencilSet, boundary::TensorGridBoundary)
+    op = boundary_restriction(g.grids[grid_id(boundary)], stencil_set, boundary_id(boundary))
+    return LazyTensors.inflate(op, size(g), grid_id(boundary))
 end
 
-"""
-    boundary_restriction(grid, stencil_set, boundary)
-
-Creates a `boundary_restriction` operator on `grid` given a `stencil_set`.
-"""
-boundary_restriction(grid, stencil_set::StencilSet, boundary) = boundary_restriction(grid, parse_stencil(stencil_set["e"]["closure"]), boundary)
+function boundary_restriction(g::EquidistantGrid, stencil_set::StencilSet, boundary)
+    closure_stencil = parse_stencil(stencil_set["e"]["closure"])
+    converted_stencil = convert(Stencil{eltype(g)}, closure_stencil)
+    return BoundaryOperator(g, converted_stencil, boundary)
+end
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,26 +1,30 @@
 """
-    normal_derivative(grid, closure_stencil::Stencil, boundary)
+    normal_derivative(g, stencil_set::StencilSet, boundary)
+    normal_derivative(g::TensorGrid, stencil_set::StencilSet, boundary::TensorGridBoundary)
+    normal_derivative(g::EquidistantGrid, stencil_set::StencilSet, boundary)
 
 Creates the normal derivative boundary operator `d` as a `LazyTensor`
 
-`d` computes the normal derivative of a grid function  on `boundary` a `Stencil` `closure_stencil`.
-`d'` is the prolongation of the normal derivative of a grid function to the whole grid using the same `closure_stencil`.
-On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
-a `BoundaryOperator`.
+`d` computes the normal derivative at `boundary` of a grid function on `g` using the
+'d1' stencil in `stencil_set`. `d'` is the prolongation of the normal
+derivative of a grid function to the whole of `g` using the same stencil. On a
+one-dimensional grid, `d` is a `BoundaryOperator`. On a multi-dimensional
+grid, `d` is the inflation of a `BoundaryOperator`.
 
 See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
-function normal_derivative(grid, closure_stencil, boundary)
-    direction = dim(boundary)
-    h_inv = inverse_spacing(grid)[direction]
+function normal_derivative end
+
 
-    op = BoundaryOperator(restrict(grid, dim(boundary)), scale(closure_stencil,h_inv), region(boundary))
-    return LazyTensors.inflate(op, size(grid), dim(boundary))
+function normal_derivative(g::TensorGrid, stencil_set::StencilSet, boundary::TensorGridBoundary)
+    op = normal_derivative(g.grids[grid_id(boundary)], stencil_set, boundary_id(boundary))
+    return LazyTensors.inflate(op, size(g), grid_id(boundary))
 end
 
-"""
-    normal_derivative(grid, stencil_set, boundary)
+function normal_derivative(g::EquidistantGrid, stencil_set::StencilSet, boundary)
+    closure_stencil = parse_stencil(stencil_set["d1"]["closure"])
+    h_inv = inverse_spacing(g)
 
-Creates a `normal_derivative` operator on `grid` given a `stencil_set`.
-"""
-normal_derivative(grid, stencil_set::StencilSet, boundary) = normal_derivative(grid, parse_stencil(stencil_set["d1"]["closure"]), boundary)
+    scaled_stencil = scale(closure_stencil,h_inv)
+    return BoundaryOperator(g, scaled_stencil, boundary)
+end
--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Sat May 20 14:33:25 2023 +0200
@@ -2,7 +2,8 @@
     ConstantInteriorScalingOperator{T,N} <: LazyTensor{T,1,1}
 
 A one-dimensional operator scaling a vector. The first and last `N` points are
-scaled with individual weights while all interior points are scaled the same.
+scaled with individual weights while all interior points are scaled using the
+same factor.
 """
 struct ConstantInteriorScalingOperator{T,N} <: LazyTensor{T,1,1}
     interior_weight::T
@@ -18,7 +19,7 @@
     end
 end
 
-function ConstantInteriorScalingOperator(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+function ConstantInteriorScalingOperator(grid::EquidistantGrid, interior_weight, closure_weights)
     return ConstantInteriorScalingOperator(interior_weight, Tuple(closure_weights), size(grid)[1])
 end
 
--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,39 +1,46 @@
 """
-    undivided_skewed04(g::EquidistantGrid, p, direction)
+    undivided_skewed04(g::TensorGrid, p, direction)
+    undivided_skewed04(g::EquidistantGrid, p)
 
-Create undivided difference operators approximating the `p`th derivative. The
-operators do not satisfy any SBP-property and are meant to be used for
+Undivided difference operators approximating the `p`th derivative. The
+operators do not satisfy any SBP property and are meant to be used for
 building artificial dissipation terms.
 
-The operators and how they are used to create accurate artifical dissipation
+The operators and how they are used to create accurate artificial dissipation
 is described in "K. Mattsson, M. Svärd, and J. Nordström, “Stable and Accurate
 Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp.
 57–79, Aug. 2004"
 """
-function undivided_skewed04(g::EquidistantGrid, p, direction)
+function undivided_skewed04 end
+
+function undivided_skewed04(g::TensorGrid, p, direction)
+    D,Dᵀ = undivided_skewed04(g.grids[direction], p)
+    return (
+        LazyTensors.inflate(D, size(g), direction),
+        LazyTensors.inflate(Dᵀ, size(g), direction),
+    )
+end
+
+function undivided_skewed04(g::EquidistantGrid, p)
     T = eltype(g)
     interior_weights = T.(dissipation_interior_weights(p))
 
-    D  = stencil_operator_distinct_closures(
+    D  = StencilOperatorDistinctClosures(
         g,
         dissipation_interior_stencil(interior_weights),
         dissipation_lower_closure_stencils(interior_weights),
         dissipation_upper_closure_stencils(interior_weights),
-        direction,
     )
-    Dᵀ = stencil_operator_distinct_closures(
+    Dᵀ = StencilOperatorDistinctClosures(
         g,
         dissipation_transpose_interior_stencil(interior_weights),
         dissipation_transpose_lower_closure_stencils(interior_weights),
         dissipation_transpose_upper_closure_stencils(interior_weights),
-        direction,
     )
 
     return D, Dᵀ
 end
 
-undivided_skewed04(g::EquidistantGrid{1}, p) = undivided_skewed04(g, p, 1)
-
 function dissipation_interior_weights(p)
    if p == 0
        return (1,)
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,48 +1,42 @@
 """
-    first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+    first_derivative(g, ..., [direction])
 
-Creates the first-derivative operator `D1` as a `LazyTensor`
+The first derivative operator `D1` as a `LazyTensor` on the given grid.
 
-`D1` approximates the first-derivative d/dξ on `grid` along the coordinate dimension specified by
-`direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions.
+`D1` approximates the first-derivative d/dξ on `g` along the coordinate
+dimension specified by `direction`.
+"""
+function first_derivative end
 
-On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the inflation of
-a `VolumeOperator`.
+"""
+    first_derivative(g::TensorGrid, stencil_set, direction)
 
 See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
-function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
-    h_inv = inverse_spacing(grid)[direction]
-
-    D₁ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd)
-    return LazyTensors.inflate(D₁, size(grid), direction)
+function first_derivative(g::TensorGrid, stencil_set, direction)
+    D₁ = first_derivative(g.grids[direction], stencil_set)
+    return LazyTensors.inflate(D₁, size(g), direction)
 end
 
-
-"""
-    first_derivative(grid, inner_stencil, closure_stencils)
-
-Creates a `first_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`.
-"""
-first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid, inner_stencil, closure_stencils, 1)
-
-
 """
-    first_derivative(grid, stencil_set::StencilSet, direction)
+    first_derivative(g::EquidistantGrid, stencil_set::StencilSet)
 
-Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`.
+The first derivative operator on an `EquidistantGrid`. 
+Uses the `D1` stencil in `stencil_set`.
 """
-function first_derivative(grid::EquidistantGrid, stencil_set::StencilSet, direction)
+function first_derivative(g::EquidistantGrid, stencil_set::StencilSet)
     inner_stencil = parse_stencil(stencil_set["D1"]["inner_stencil"])
     closure_stencils = parse_stencil.(stencil_set["D1"]["closure_stencils"])
-    first_derivative(grid,inner_stencil,closure_stencils,direction);
+    return first_derivative(g, inner_stencil, closure_stencils);
 end
 
-
 """
-    first_derivative(grid, stencil_set)
+    first_derivative(g::EquidistantGrid, inner_stencil::Stencil, closure_stencils)
 
-Creates a `first_derivative` operator on a 1D `grid` given a `stencil_set`.
+The first derivative operator on an `EquidistantGrid` given an
+`inner_stencil` and `closure_stencils`.
 """
-first_derivative(grid::EquidistantGrid{1}, stencil_set::StencilSet) = first_derivative(grid, stencil_set, 1)
+function first_derivative(g::EquidistantGrid, inner_stencil::Stencil, closure_stencils)
+    h⁻¹ = inverse_spacing(g)
+    return VolumeOperator(g, scale(inner_stencil,h⁻¹), scale.(closure_stencils,h⁻¹), odd)
+end
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,48 +1,37 @@
 """
-    second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
-
-Creates the second-derivative operator `D2` as a `LazyTensor`
+    second_derivative(g::EquidistantGrid, stencil_set, direction)
 
-`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by
-`direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions.
+Creates the second derivative operator `D2` as a `LazyTensor`
 
-On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the inflation of
-a `VolumeOperator`.
+`D2` approximates the second-derivative d²/dξ² on `g` along the coordinate
+dimension specified by `direction`.
 
 See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
-function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
-    h_inv = inverse_spacing(grid)[direction]
-
-    D₂ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even)
-    return LazyTensors.inflate(D₂, size(grid), direction)
+function second_derivative(g::TensorGrid, stencil_set, direction)
+    D₂ = second_derivative(g.grids[direction], stencil_set)
+    return LazyTensors.inflate(D₂, size(g), direction)
 end
 
-
-"""
-    second_derivative(grid, inner_stencil, closure_stencils)
-
-Creates a `second_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`.
-"""
-second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid, inner_stencil, closure_stencils,1)
-
-
 """
-    second_derivative(grid, stencil_set, direction)
+    second_derivative(g::EquidistantGrid, stencil_set::::StencilSet)
 
-Creates a `second_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`.
+The second derivative operator on an `EquidistantGrid`. 
+Uses the `D2` stencil in `stencil_set`.
 """
-function second_derivative(grid::EquidistantGrid, stencil_set::StencilSet, direction)
+function second_derivative(g::EquidistantGrid, stencil_set::StencilSet)
     inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
     closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-    second_derivative(grid,inner_stencil,closure_stencils,direction);
+    return second_derivative(g, inner_stencil, closure_stencils)
 end
 
-
 """
-    second_derivative(grid, stencil_set)
+    second_derivative(g::EquidistantGrid, inner_stencil::Stencil, closure_stencils)
 
-Creates a `second_derivative` operator on a 1D `grid` given a `stencil_set`.
+The second derivative operator on an `EquidistantGrid`, given `inner_stencil` and
+`closure_stencils`.
 """
-second_derivative(grid::EquidistantGrid{1}, stencil_set::StencilSet) = second_derivative(grid, stencil_set, 1)
\ No newline at end of file
+function second_derivative(g::EquidistantGrid, inner_stencil::Stencil, closure_stencils)
+    h⁻¹ = inverse_spacing(g)
+    return VolumeOperator(g, scale(inner_stencil,h⁻¹^2), scale.(closure_stencils,h⁻¹^2), even)
+end
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,46 +1,52 @@
 """
-    inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
+    inner_product(grid, ...)
 
-Creates the discrete inner product operator `H` as a `LazyTensor` on an
-equidistant grid, defined as `(u,v)  = u'Hv` for grid functions `u,v`.
+The inner product on a given grid with weights from a stencils set or given
+explicitly.
+"""
+function inner_product end
+
+"""
+    inner_product(tg::TensorGrid, stencil_set::StencilSet)
 
-`inner_product` creates `H` on `grid` using the `interior_weight` for the
-interior points and the `closure_weights` for the points close to the
-boundary.
+The inner product on `tg`, i.e., the tensor product of the
+individual grids' inner products, using weights `H` from `stencil_set`.
+"""
+function inner_product(tg::TensorGrid, stencil_set::StencilSet)
+    return mapreduce(g->inner_product(g,stencil_set), ⊗, tg.grids)
+end
 
-On a 1-dimensional grid, `H` is a `ConstantInteriorScalingOperator`. On a
-N-dimensional grid, `H` is the outer product of the 1-dimensional inner
-product operators for each coordinate direction. On a 0-dimensional grid,
-`H` is a 0-dimensional `IdentityTensor`.
+"""
+    inner_product(g::EquidistantGrid, stencil_set::StencilSet)
+
+The inner product on `g` using weights `H` from `stencil_set`.
 
 See also: [`ConstantInteriorScalingOperator`](@ref).
 """
-function inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
-    Hs = ()
-
-    for i ∈ dims(grid)
-        Hs = (Hs..., inner_product(restrict(grid, i), interior_weight, closure_weights))
-    end
-
-    return foldl(⊗, Hs)
+function inner_product(g::EquidistantGrid, stencil_set::StencilSet)
+    interior_weight = parse_scalar(stencil_set["H"]["inner"])
+    closure_weights = parse_tuple(stencil_set["H"]["closure"])
+    return inner_product(g, interior_weight, closure_weights)
 end
 
-function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
-    h = spacing(grid)[1]
+"""
+    inner_product(g::EquidistantGrid, interior_weight, closure_weights)
+
+The inner product on `g` with explicit weights.
 
-    H = SbpOperators.ConstantInteriorScalingOperator(grid, h*interior_weight, h.*closure_weights)
-    return H
+See also: [`ConstantInteriorScalingOperator`](@ref).
+"""
+function inner_product(g::EquidistantGrid, interior_weight, closure_weights)
+    h = spacing(g)
+    return SbpOperators.ConstantInteriorScalingOperator(g, h*interior_weight, h.*closure_weights)
 end
 
-inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityTensor{eltype(grid)}()
-
 """
-    inner_product(grid, stencil_set)
+    inner_product(g::ZeroDimGrid, stencil_set::StencilSet)
 
-Creates a `inner_product` operator on `grid` given a `stencil_set`.
+The identity tensor with the correct type parameters.
+
+Implemented to simplify 1D code for SBP operators.
 """
-function inner_product(grid, stencil_set::StencilSet)
-    inner_stencil = parse_scalar(stencil_set["H"]["inner"])
-    closure_stencils = parse_tuple(stencil_set["H"]["closure"])
-    return inner_product(grid, inner_stencil, closure_stencils)
-end
+inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}()
+
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,42 +1,51 @@
 """
-    inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
+    inverse_inner_product(grid, ...)
+
+The inverse inner product on a given grid with weights from a stencils set or given
+explicitly.
+"""
+function inverse_inner_product end
+
+"""
+    inverse_inner_product(tg::TensorGrid, stencil_set::StencilSet)
 
-Constructs the inverse inner product operator `H⁻¹` as a `LazyTensor` using
-the weights of `H`, `interior_weight`, `closure_weights`. `H⁻¹` is inverse of
-the inner product operator `H`.
+The inverse of inner product on `tg`, i.e., the tensor product of the
+individual grids' inverse inner products, using weights `H` from `stencil_set`.
+"""
+function inverse_inner_product(tg::TensorGrid, stencil_set::StencilSet)
+    return mapreduce(g->inverse_inner_product(g,stencil_set), ⊗, tg.grids)
+end
 
-On a 1-dimensional grid, `H⁻¹` is a `ConstantInteriorScalingOperator`. On an
-N-dimensional grid, `H⁻¹` is the outer product of the 1-dimensional inverse
-inner product operators for each coordinate direction. On a 0-dimensional
-`grid`, `H⁻¹` is a 0-dimensional `IdentityTensor`.
+"""
+    inverse_inner_product(g::EquidistantGrid, stencil_set::StencilSet)
+
+The inverse of the inner product on `g` using weights `H` from `stencil_set`.
 
 See also: [`ConstantInteriorScalingOperator`](@ref).
 """
-function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
-    H⁻¹s = ()
-
-    for i ∈ dims(grid)
-        H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), interior_weight, closure_weights))
-    end
-
-    return foldl(⊗, H⁻¹s)
+function inverse_inner_product(g::EquidistantGrid, stencil_set::StencilSet)
+    interior_weight = parse_scalar(stencil_set["H"]["inner"])
+    closure_weights = parse_tuple(stencil_set["H"]["closure"])
+    return inverse_inner_product(g, interior_weight, closure_weights)
 end
 
-function inverse_inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
-    h⁻¹ = inverse_spacing(grid)[1]
-    H⁻¹ = SbpOperators.ConstantInteriorScalingOperator(grid, h⁻¹*1/interior_weight, h⁻¹./closure_weights)
-    return H⁻¹
+"""
+    inverse_inner_product(g::EquidistantGrid, interior_weight, closure_weights)
+
+The inverse inner product on `g` with explicit weights.
+
+See also: [`ConstantInteriorScalingOperator`](@ref).
+"""
+function inverse_inner_product(g::EquidistantGrid, interior_weight, closure_weights)
+    h⁻¹ = inverse_spacing(g)
+    return SbpOperators.ConstantInteriorScalingOperator(g, h⁻¹*1/interior_weight, h⁻¹./closure_weights)
 end
 
-inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityTensor{eltype(grid)}()
-
 """
-    inverse_inner_product(grid, stencil_set)
+    inverse_inner_product(g::ZeroDimGrid, stencil_set::StencilSet)
 
-Creates a `inverse_inner_product` operator on `grid` given a `stencil_set`.
+The identity tensor with the correct type parameters.
+
+Implemented to simplify 1D code for SBP operators.
 """
-function inverse_inner_product(grid, stencil_set::StencilSet)
-    inner_stencil = parse_scalar(stencil_set["H"]["inner"])
-    closure_stencils = parse_tuple(stencil_set["H"]["closure"])
-    return inverse_inner_product(grid, inner_stencil, closure_stencils)
-end
+inverse_inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}()
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,9 +1,8 @@
 """
     Laplace{T, Dim, TM} <: LazyTensor{T, Dim, Dim}
 
-Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
-`LazyTensor`. Additionally `Laplace` stores the `StencilSet`
-used to construct the `LazyTensor `.
+The Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
+`LazyTensor`.
 """
 struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim}
     D::TM       # Difference operator
@@ -11,17 +10,15 @@
 end
 
 """
-    Laplace(grid::Equidistant, stencil_set)
+    Laplace(g::Grid, stencil_set::StencilSet)
 
-Creates the `Laplace` operator `Δ` on `grid` given a `stencil_set`. 
+Creates the `Laplace` operator `Δ` on `g` given `stencil_set`. 
 
 See also [`laplace`](@ref).
 """
-function Laplace(grid::EquidistantGrid, stencil_set::StencilSet)
-    inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-    Δ = laplace(grid, inner_stencil,closure_stencils)
-    return Laplace(Δ,stencil_set)
+function Laplace(g::Grid, stencil_set::StencilSet)
+    Δ = laplace(g, stencil_set)
+    return Laplace(Δ, stencil_set)
 end
 
 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
@@ -32,24 +29,26 @@
 # Base.show(io::IO, L::Laplace) = ...
 
 """
-    laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
-
-Creates the Laplace operator operator `Δ` as a `LazyTensor`
+    laplace(g::Grid, stencil_set)
 
-`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using
-the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions.
+Creates the Laplace operator operator `Δ` as a `LazyTensor` on `g`.
 
-On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
-multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
-where the sum is carried out lazily.
+`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `g`. The
+approximation depends on the type of grid and the stencil set.
 
 See also: [`second_derivative`](@ref).
 """
-function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
-    Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
-    for d = 2:ndims(grid)
-        Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
+function laplace end
+function laplace(g::TensorGrid, stencil_set)
+    # return mapreduce(+, enumerate(g.grids)) do (i, gᵢ)
+    #     Δᵢ = laplace(gᵢ, stencil_set)
+    #     LazyTensors.inflate(Δᵢ, size(g), i)
+    # end
+
+    Δ = LazyTensors.inflate(laplace(g.grids[1], stencil_set), size(g), 1)
+    for d = 2:ndims(g)
+        Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d)
     end
     return Δ
 end
+laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
--- a/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,26 +1,8 @@
-"""
-    stencil_operator_distinct_closures(
-        grid::EquidistantGrid,
-        inner_stencil,
-        lower_closure,
-        upper_closure,
-        direction
-    )
-
-Creates a multi-dimensional `StencilOperatorDistinctClosures` acting on grid
-functions of `grid`.
-
-See also: [`StencilOperatorDistinctClosures`](@ref)
-"""
-function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction)
-    op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure)
-    return LazyTensors.inflate(op, size(grid), direction)
-end
-
 """
     StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1}
 
-A one dimensional stencil operator with separate closures for the two boundaries.
+A one dimensional stencil operator with separate closures for the two
+boundaries.
 
 `StencilOperatorDistinctClosures` can be contrasted to `VolumeOperator` in
 that it has different closure stencils for the upper and lower boundary.
@@ -28,7 +10,7 @@
 closures is useful for representing operators with skewed stencils like upwind
 operators.
 
-See also: [`VolumeOperator`](@ref), [`stencil_operator_distinct_closures`](@ref)
+See also: [`VolumeOperator`](@ref)
 """
 struct StencilOperatorDistinctClosures{T,K,N,M,LC<:NTuple{N,Stencil{T,L}} where L, UC<:NTuple{M,Stencil{T,L}} where L} <: LazyTensor{T,1,1}
     inner_stencil::Stencil{T,K}
@@ -37,7 +19,7 @@
     size::Tuple{Int}
 end
 
-function StencilOperatorDistinctClosures(grid::EquidistantGrid{1}, inner_stencil, lower_closure, upper_closure)
+function StencilOperatorDistinctClosures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure)
     return StencilOperatorDistinctClosures(inner_stencil, Tuple(lower_closure), Tuple(upper_closure), size(grid))
 end
 
--- a/src/SbpOperators/volumeops/volume_operator.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,7 +1,7 @@
 """
     VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
 
-Implements a one-dimensional constant coefficients volume operator
+A one-dimensional constant coefficients stencil operator.
 """
 struct VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
     inner_stencil::Stencil{T,N}
@@ -10,7 +10,7 @@
     parity::Parity
 end
 
-function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity)
+function VolumeOperator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity)
     return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity)
 end
 
--- a/test/Grids/equidistant_grid_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -1,39 +1,114 @@
 using Sbplib.Grids
 using Test
 using Sbplib.RegionIndices
+using Sbplib.LazyTensors
 
 
 @testset "EquidistantGrid" begin
-    @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid
-    @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid
-    # constuctor
-    @test_throws DomainError EquidistantGrid(0,0.0,1.0)
-    @test_throws DomainError EquidistantGrid(1,1.0,1.0)
-    @test_throws DomainError EquidistantGrid(1,1.0,-1.0)
-    @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,))
+    @test EquidistantGrid(0:0.1:10) isa EquidistantGrid
+    @test EquidistantGrid(range(0,1,length=10)) isa EquidistantGrid
+    @test EquidistantGrid(LinRange(0,1,11)) isa EquidistantGrid
+
+    @testset "Indexing Interface" begin
+        g = EquidistantGrid(0:0.1:10)
+        @test g[1] == 0.0
+        @test g[5] == 0.4
+        @test g[101] == 10.0
+
+        @test g[begin] == 0.0
+        @test g[end] == 10.0
+
+        @test all(eachindex(g) .== 1:101)
+    end
+
+    @testset "Iterator interface" begin
+        @test eltype(EquidistantGrid(0:10)) == Int
+        @test eltype(EquidistantGrid(0:2:10)) == Int
+        @test eltype(EquidistantGrid(0:0.1:10)) == Float64
+        @test size(EquidistantGrid(0:10)) == (11,)
+        @test size(EquidistantGrid(0:0.1:10)) == (101,)
+
+        @test collect(EquidistantGrid(0:0.1:0.5)) == [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]
+
+        @test Base.IteratorSize(EquidistantGrid{Float64, StepRange{Float64}}) == Base.HasShape{1}()
+    end
 
     @testset "Base" begin
-        @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64
-        @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int
-        @test size(EquidistantGrid(4,0.0,1.0)) == (4,)
-        @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
-        @test ndims(EquidistantGrid(4,0.0,1.0)) == 1
-        @test ndims(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == 2
+        @test ndims(EquidistantGrid(0:10)) == 1
     end
 
     @testset "spacing" begin
-        @test [spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(1. /3,)...] atol=5e-13
-        @test [spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(0.5, 1.)...] atol=5e-13
+        @test spacing(EquidistantGrid(0:10)) == 1
+        @test spacing(EquidistantGrid(0:0.1:10)) == 0.1
     end
 
     @testset "inverse_spacing" begin
-        @test [inverse_spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(3.,)...] atol=5e-13
-        @test [inverse_spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(2, 1.)...] atol=5e-13
+        @test inverse_spacing(EquidistantGrid(0:10)) == 1
+        @test inverse_spacing(EquidistantGrid(0:0.1:10)) == 10
+    end
+
+    @testset "boundary_identifiers" begin
+        g = EquidistantGrid(0:0.1:10)
+        @test boundary_identifiers(g) == (Lower(), Upper())
+        @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)
+    end
+
+    @testset "refine" begin
+        g = EquidistantGrid(0:0.1:1)
+        @test refine(g, 1) == g
+        @test refine(g, 2) == EquidistantGrid(0:0.05:1)
+        @test refine(g, 3) == EquidistantGrid(0:(0.1/3):1)
     end
 
-    @testset "points" begin
-        g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11))
-        gp = points(g);
+    @testset "coarsen" begin
+        g = EquidistantGrid(0:1:10)
+        @test coarsen(g, 1) == g
+        @test coarsen(g, 2) == EquidistantGrid(0:2:10)
+
+        g = EquidistantGrid(0:0.1:1)
+        @test coarsen(g, 1) == g
+        @test coarsen(g, 2) == EquidistantGrid(0:0.2:1)
+
+        g = EquidistantGrid(0:10)
+        @test coarsen(g, 1) == EquidistantGrid(0:1:10)
+        @test coarsen(g, 2) == EquidistantGrid(0:2:10)
+
+        @test_throws DomainError(3, "Size minus 1 must be divisible by the ratio.") coarsen(g, 3)
+    end
+end
+
+
+@testset "equidistant_grid" begin
+    @test equidistant_grid(4,0.0,1.0) isa EquidistantGrid
+    @test equidistant_grid((4,3),(0.0,0.0),(8.0,5.0)) isa TensorGrid
+
+    # constuctor
+    @test_throws DomainError equidistant_grid(0,0.0,1.0)
+    @test_throws DomainError equidistant_grid(1,1.0,1.0)
+    @test_throws DomainError equidistant_grid(1,1.0,-1.0)
+
+    @test_throws DomainError equidistant_grid((0,0),(0.0,0.0),(1.0,1.0))
+    @test_throws DomainError equidistant_grid((1,1),(1.0,1.0),(1.0,1.0))
+    @test_throws DomainError equidistant_grid((1,1),(1.0,1.0),(-1.0,-1.0))
+
+    @testset "Base" begin
+        @test eltype(equidistant_grid(4,0.0,1.0)) == Float64
+        @test eltype(equidistant_grid((4,3),(0,0),(1,3))) <: AbstractVector{Float64}
+        @test size(equidistant_grid(4,0.0,1.0)) == (4,)
+        @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+        @test ndims(equidistant_grid(4,0.0,1.0)) == 1
+        @test ndims(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == 2
+    end
+
+    @testset "getindex" begin
+        g = equidistant_grid((5,3), (-1.0,0.0), (0.0,7.11))
+        gp = collect(g);
         p = [(-1.,0.)      (-1.,7.11/2)   (-1.,7.11);
             (-0.75,0.)    (-0.75,7.11/2) (-0.75,7.11);
             (-0.5,0.)     (-0.5,7.11/2)  (-0.5,7.11);
@@ -43,83 +118,18 @@
             @test [gp[i]...] ≈ [p[i]...] atol=5e-13
         end
     end
+end
 
-    @testset "restrict" begin
-        g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))
-        @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0)
-        @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0)
 
-        g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
-        @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0)
-        @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0)
-        @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0)
-        @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0))
-        @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
-        @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0))
-        @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0))
-    end
-
-    @testset "boundary_identifiers" begin
-        g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
-        bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(),
-                CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(),
-                CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}())
-        @test boundary_identifiers(g) == bids
-        @inferred boundary_identifiers(g)
-    end
+@testset "change_length" begin
+    @test Grids.change_length(0:20, 21) == 0:20
+    @test Grids.change_length(0:20, 11) == 0:2:20
+    @test Grids.change_length(0:2:20, 21) == 0:20
 
-    @testset "boundary_grid" begin
-            @testset "1D" begin
-                g = EquidistantGrid(5,0.0,2.0)
-                (id_l, id_r) = boundary_identifiers(g)
-                @test boundary_grid(g,id_l) == EquidistantGrid{Float64}()
-                @test boundary_grid(g,id_r) == EquidistantGrid{Float64}()
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}())
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}())
-            end
-            @testset "2D" begin
-                g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
-                (id_w, id_e, id_s, id_n) = boundary_identifiers(g)
-                @test boundary_grid(g,id_w) == restrict(g,2)
-                @test boundary_grid(g,id_e) == restrict(g,2)
-                @test boundary_grid(g,id_s) == restrict(g,1)
-                @test boundary_grid(g,id_n) == restrict(g,1)
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
-            end
-            @testset "3D" begin
-                g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
-                (id_w, id_e,
-                 id_s, id_n,
-                 id_t, id_b) = boundary_identifiers(g)
-                @test boundary_grid(g,id_w) == restrict(g,[2,3])
-                @test boundary_grid(g,id_e) == restrict(g,[2,3])
-                @test boundary_grid(g,id_s) == restrict(g,[1,3])
-                @test boundary_grid(g,id_n) == restrict(g,[1,3])
-                @test boundary_grid(g,id_t) == restrict(g,[1,2])
-                @test boundary_grid(g,id_b) == restrict(g,[1,2])
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
-            end
-    end
+    @test Grids.change_length(range(0,1,length=10), 10) == range(0,1,length=10)
+    @test Grids.change_length(range(0,1,length=10), 5) == range(0,1,length=5)
+    @test Grids.change_length(range(0,1,length=10), 20) == range(0,1,length=20)
 
-    @testset "refine" begin
-        @test refine(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}()
-        @test refine(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}()
-
-        g = EquidistantGrid((10,5),(0.,1.),(2.,3.))
-        @test refine(g, 1) == g
-        @test refine(g, 2) == EquidistantGrid((19,9),(0.,1.),(2.,3.))
-        @test refine(g, 3) == EquidistantGrid((28,13),(0.,1.),(2.,3.))
-    end
-
-    @testset "coarsen" begin
-        @test coarsen(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}()
-        @test coarsen(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}()
-
-        g = EquidistantGrid((7,13),(0.,1.),(2.,3.))
-        @test coarsen(g, 1) == g
-        @test coarsen(g, 2) == EquidistantGrid((4,7),(0.,1.),(2.,3.))
-        @test coarsen(g, 3) == EquidistantGrid((3,5),(0.,1.),(2.,3.))
-
-        @test_throws DomainError(4, "Size minus 1 must be divisible by the ratio.") coarsen(g, 4) == EquidistantGrid((3,5),(0.,1.),(2.,3.))
-    end
+    @test Grids.change_length(LinRange(1,2,10),10) == LinRange(1,2,10)
+    @test Grids.change_length(LinRange(1,2,10),15) == LinRange(1,2,15)
 end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Grids/grid_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -0,0 +1,69 @@
+using Test
+using Sbplib.Grids
+using Sbplib.LazyTensors
+using StaticArrays
+
+@testset "Grid" begin
+    struct DummyGrid{T,D} <: Grid{T,D} end
+
+    @test eltype(DummyGrid{Int, 2}) == Int
+    @test eltype(DummyGrid{Int, 2}()) == Int
+
+    @test ndims(DummyGrid{Int, 2}()) == 2
+
+    @test coordinate_size(DummyGrid{Int, 1}()) == 1
+    @test coordinate_size(DummyGrid{SVector{3,Float64}, 2}()) == 3
+
+    @test coordinate_size(DummyGrid{SVector{3,Float64}, 2}) == 3
+
+    @testset "component_type" begin
+        @test component_type(DummyGrid{Int,1}()) == Int
+        @test component_type(DummyGrid{Float64,1}()) == Float64
+        @test component_type(DummyGrid{Rational,1}()) == Rational
+
+        @test component_type(DummyGrid{SVector{3,Int},2}()) == Int
+        @test component_type(DummyGrid{SVector{2,Float64},3}()) == Float64
+        @test component_type(DummyGrid{SVector{4,Rational},4}()) == Rational
+
+        @test component_type(DummyGrid{Float64,1}) == Float64
+        @test component_type(DummyGrid{SVector{2,Float64},3}) == Float64
+    end
+end
+
+@testset "eval_on" begin
+    @test eval_on(ZeroDimGrid(@SVector[1.,2.]), x̄->x̄[1]+x̄[2]) isa LazyArray
+    @test eval_on(ZeroDimGrid(@SVector[1.,2.]), x̄->x̄[1]+x̄[2]) == fill(3.)
+    @test eval_on(ZeroDimGrid(@SVector[3.,2.]), x̄->x̄[1]+x̄[2]) == fill(5.)
+
+    @test eval_on(ZeroDimGrid(1.), x̄->2x̄) isa LazyArray
+    @test eval_on(ZeroDimGrid(1.), x̄->2x̄) == fill(2.)
+
+    @test eval_on(EquidistantGrid(range(0,1,length=4)), x->2x) isa LazyArray
+    @test eval_on(EquidistantGrid(range(0,1,length=4)), x->2x) == 2 .* range(0,1,length=4)
+
+
+    g = equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))
+
+    @test eval_on(g, x̄ -> 0.) isa LazyArray
+    @test eval_on(g, x̄ -> 0.) == fill(0., (5,3))
+
+    @test eval_on(g, x̄ -> sin(x̄[1])*cos(x̄[2])) == map(x̄->sin(x̄[1])*cos(x̄[2]), g)
+
+    # Vector valued function
+    @test eval_on(g, x̄ -> @SVector[x̄[2], x̄[1]]) isa LazyArray{SVector{2,Float64}}
+    @test eval_on(g, x̄ -> @SVector[x̄[2], x̄[1]]) == map(x̄ -> @SVector[x̄[2], x̄[1]], g)
+
+    # Multi-argument functions
+    f(x,y) = sin(x)*cos(y)
+    @test eval_on(g, f) == map(x̄->f(x̄...), g)
+end
+
+@testset "_ncomponents" begin
+    @test Grids._ncomponents(Int) == 1
+    @test Grids._ncomponents(Float64) == 1
+    @test Grids._ncomponents(Rational) == 1
+
+    @test Grids._ncomponents(SVector{3,Int}) == 3
+    @test Grids._ncomponents(SVector{2,Float64}) == 2
+    @test Grids._ncomponents(SVector{4,Rational}) == 4
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Grids/tensor_grid_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -0,0 +1,146 @@
+using Test
+using Sbplib.Grids
+using StaticArrays
+using Sbplib.RegionIndices
+
+@testset "TensorGrid" begin
+    g₁ = EquidistantGrid(range(0,1,length=11))
+    g₂ = EquidistantGrid(range(2,3,length=6))
+    g₃ = EquidistantGrid(1:10)
+    g₄ = ZeroDimGrid(@SVector[1,2])
+
+    @test TensorGrid(g₁, g₂) isa TensorGrid
+    @test TensorGrid(g₁, g₂) isa Grid{SVector{2,Float64}, 2}
+    @test TensorGrid(g₃, g₃) isa Grid{SVector{2,Int}, 2}
+    @test TensorGrid(g₁, g₂, g₃) isa Grid{SVector{3,Float64}, 3}
+    @test TensorGrid(g₁, g₄) isa Grid{SVector{3,Float64}, 1}
+    @test TensorGrid(g₁, g₄, g₂) isa Grid{SVector{4,Float64}, 2}
+
+    @testset "Indexing Interface" begin
+        @testset "regular indexing" begin
+            @test TensorGrid(g₁, g₂)[1,1] isa SVector{2,Float64}
+            @test TensorGrid(g₁, g₂)[1,1] == [0.0,2.0]
+            @test TensorGrid(g₁, g₂)[3,5] == [0.2,2.8]
+            @test TensorGrid(g₁, g₂)[10,6] == [0.9,3.0]
+
+            @test TensorGrid(g₁, g₃)[1,1] isa SVector{2,Float64}
+            @test TensorGrid(g₁, g₃)[1,1] == [0.0,1.0]
+
+            @test TensorGrid(g₁, g₂, g₃)[3,4,5] isa SVector{3,Float64}
+            @test TensorGrid(g₁, g₂, g₃)[3,4,5] == [0.2, 2.6, 5.0]
+
+            @test TensorGrid(g₁, g₄)[3] isa SVector{3,Float64}
+            @test TensorGrid(g₁, g₄)[3] == [0.2, 1., 2.]
+
+            @test TensorGrid(g₁, g₄, g₂)[3,2] isa SVector{4,Float64}
+            @test TensorGrid(g₁, g₄, g₂)[3,2] == [0.2, 1., 2., 2.2]
+        end
+
+        @testset "cartesian indexing" begin
+            cases = [
+                (TensorGrid(g₁, g₂),     (1,1)  ),
+                (TensorGrid(g₁, g₂),     (3,5)  ),
+                (TensorGrid(g₁, g₂),     (10,6) ),
+                (TensorGrid(g₁, g₃),     (1,1)  ),
+                (TensorGrid(g₁, g₂, g₃), (3,4,5)),
+                (TensorGrid(g₁, g₄),     (3)    ),
+                (TensorGrid(g₁, g₄, g₂), (3,2)  ),
+            ]
+
+            @testset "i = $is" for (g, is) ∈ cases
+                @test g[CartesianIndex(is...)] == g[is...]
+            end
+        end
+
+        @testset "eachindex" begin
+            @test eachindex(TensorGrid(g₁, g₂)) == CartesianIndices((11,6))
+            @test eachindex(TensorGrid(g₁, g₃)) == CartesianIndices((11,10))
+            @test eachindex(TensorGrid(g₁, g₂, g₃)) == CartesianIndices((11,6,10))
+            @test eachindex(TensorGrid(g₁, g₄)) == CartesianIndices((11,))
+            @test eachindex(TensorGrid(g₁, g₄, g₂)) == CartesianIndices((11,6))
+        end
+    end
+
+    @testset "Iterator interface" begin
+        @test eltype(TensorGrid(g₁, g₂)) == SVector{2,Float64}
+        @test eltype(TensorGrid(g₁, g₃)) == SVector{2,Float64}
+        @test eltype(TensorGrid(g₁, g₂, g₃)) == SVector{3,Float64}
+        @test eltype(TensorGrid(g₁, g₄)) == SVector{3,Float64}
+        @test eltype(TensorGrid(g₁, g₄, g₂)) == SVector{4,Float64}
+
+        @test size(TensorGrid(g₁, g₂)) == (11,6)
+        @test size(TensorGrid(g₁, g₃)) == (11,10)
+        @test size(TensorGrid(g₁, g₂, g₃)) == (11,6,10)
+        @test size(TensorGrid(g₁, g₄)) == (11,)
+        @test size(TensorGrid(g₁, g₄, g₂)) == (11,6)
+
+        @test Base.IteratorSize(TensorGrid(g₁, g₂)) == Base.HasShape{2}()
+        @test Base.IteratorSize(TensorGrid(g₁, g₃)) == Base.HasShape{2}()
+        @test Base.IteratorSize(TensorGrid(g₁, g₂, g₃)) == Base.HasShape{3}()
+        @test Base.IteratorSize(TensorGrid(g₁, g₄)) == Base.HasShape{1}()
+        @test Base.IteratorSize(TensorGrid(g₁, g₄, g₂)) == Base.HasShape{2}()
+
+        @test iterate(TensorGrid(g₁, g₂))[1] isa SVector{2,Float64}
+        @test iterate(TensorGrid(g₁, g₃))[1] isa SVector{2,Float64}
+        @test iterate(TensorGrid(g₁, g₂, g₃))[1] isa SVector{3,Float64}
+        @test iterate(TensorGrid(g₁, g₄))[1] isa SVector{3,Float64}
+        @test iterate(TensorGrid(g₁, g₄, g₂))[1] isa SVector{4,Float64}
+
+        @test collect(TensorGrid(g₁, g₂)) == [@SVector[x,y] for x ∈ range(0,1,length=11), y ∈ range(2,3,length=6)]
+        @test collect(TensorGrid(g₁, g₃)) == [@SVector[x,y] for x ∈ range(0,1,length=11), y ∈ 1:10]
+        @test collect(TensorGrid(g₁, g₂, g₃)) == [@SVector[x,y,z] for x ∈ range(0,1,length=11), y ∈ range(2,3,length=6), z ∈ 1:10]
+        @test collect(TensorGrid(g₁, g₄)) == [@SVector[x,1,2] for x ∈ range(0,1,length=11)]
+        @test collect(TensorGrid(g₁, g₄, g₂)) == [@SVector[x,1,2,y] for x ∈ range(0,1,length=11), y ∈ range(2,3,length=6)]
+    end
+
+    @testset "refine" begin
+        g1(n) = EquidistantGrid(range(0,1,length=n))
+        g2(n) = EquidistantGrid(range(2,3,length=n))
+
+        @test refine(TensorGrid(g1(11), g2(6)),1) == TensorGrid(g1(11), g2(6))
+        @test refine(TensorGrid(g1(11), g2(6)),2) == TensorGrid(g1(21), g2(11))
+        @test refine(TensorGrid(g1(11), g2(6)),3) == TensorGrid(g1(31), g2(16))
+        @test refine(TensorGrid(g1(11), g₄), 1) == TensorGrid(g1(11), g₄)
+        @test refine(TensorGrid(g1(11), g₄), 2) == TensorGrid(g1(21), g₄)
+    end
+
+    @testset "coarsen" begin
+        g1(n) = EquidistantGrid(range(0,1,length=n))
+        g2(n) = EquidistantGrid(range(2,3,length=n))
+
+        @test coarsen(TensorGrid(g1(11), g2(6)),1) == TensorGrid(g1(11), g2(6))
+        @test coarsen(TensorGrid(g1(21), g2(11)),2) == TensorGrid(g1(11), g2(6))
+        @test coarsen(TensorGrid(g1(31), g2(16)),3) == TensorGrid(g1(11), g2(6))
+        @test coarsen(TensorGrid(g1(11), g₄), 1) == TensorGrid(g1(11), g₄)
+        @test coarsen(TensorGrid(g1(21), g₄), 2) == TensorGrid(g1(11), g₄)
+    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}())
+    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₄)
+    end
+end
+
+@testset "combined_coordinate_vector_type" begin
+    @test Grids.combined_coordinate_vector_type(Float64) == Float64
+    @test Grids.combined_coordinate_vector_type(Float64, Int) == SVector{2,Float64}
+    @test Grids.combined_coordinate_vector_type(Float32, Int16, Int32) == SVector{3,Float32}
+
+    @test Grids.combined_coordinate_vector_type(SVector{2,Float64}) == SVector{2,Float64}
+    @test Grids.combined_coordinate_vector_type(SVector{2,Float64}, SVector{1,Float64}) == SVector{3,Float64}
+    @test Grids.combined_coordinate_vector_type(SVector{2,Float64}, SVector{1,Int}, SVector{3, Float32}) == SVector{6,Float64}
+end
+
+@testset "combine_coordinates" begin
+    @test Grids.combine_coordinates(1,2,3) isa SVector{3, Int}
+    @test Grids.combine_coordinates(1,2,3) == [1,2,3]
+    @test Grids.combine_coordinates(1,2.,3) isa SVector{3, Float64}
+    @test Grids.combine_coordinates(1,2.,3) == [1,2,3]
+    @test Grids.combine_coordinates(1,@SVector[2.,3]) isa SVector{3, Float64}
+    @test Grids.combine_coordinates(1,@SVector[2.,3]) == [1,2,3]
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Grids/zero_dim_grid_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -0,0 +1,44 @@
+using Test
+using Sbplib.Grids
+using StaticArrays
+
+@testset "ZeroDimGrid" begin
+    @test ZeroDimGrid(1) isa ZeroDimGrid{Int}
+    @test ZeroDimGrid([1,2,3]) isa ZeroDimGrid{Vector{Int}}
+    @test ZeroDimGrid(@SVector[1.0,2.0]) isa ZeroDimGrid{SVector{2,Float64}}
+
+    @testset "Indexing Interface" begin
+        g = ZeroDimGrid(@SVector[1,2])
+
+        @test g[] == [1,2]
+        @test eachindex(g) == CartesianIndices(())
+    end
+
+    @testset "Iterator interface" begin
+        g = ZeroDimGrid(@SVector[1,2])
+
+        @test Base.IteratorSize(g) == Base.HasShape{0}()
+        @test eltype(g) == SVector{2,Int}
+        @test length(g) == 1
+        @test size(g) == ()
+        @test collect(g) == fill(@SVector[1,2])
+    end
+
+    @testset "refine" begin
+        @test refine(ZeroDimGrid(@SVector[1.0,2.0]),1) == ZeroDimGrid(@SVector[1.0,2.0])
+        @test refine(ZeroDimGrid(@SVector[1.0,2.0]),2) == ZeroDimGrid(@SVector[1.0,2.0])
+    end
+
+    @testset "coarsen" begin
+        @test coarsen(ZeroDimGrid(@SVector[1.0,2.0]),1) == ZeroDimGrid(@SVector[1.0,2.0])
+        @test coarsen(ZeroDimGrid(@SVector[1.0,2.0]),2) == ZeroDimGrid(@SVector[1.0,2.0])
+    end
+
+    @testset "boundary_identifiers" begin
+        @test boundary_identifiers(ZeroDimGrid(@SVector[1.0,2.0])) == ()
+    end
+
+    @testset "boundary_grid" begin
+        @test_throws ArgumentError("ZeroDimGrid has no boundaries") boundary_grid(ZeroDimGrid(@SVector[1.0,2.0]), :bid)
+    end
+end
--- a/test/Manifest.toml	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/Manifest.toml	Sat May 20 14:33:25 2023 +0200
@@ -2,7 +2,7 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "23260eda65ade7d11fffed313a68520d0bc053fc"
+project_hash = "f2b0634c12bbed93a17efc88d466604d5a07c465"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -250,6 +250,17 @@
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.1.7"
 
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
+git-tree-sha1 = "2d7d9e1ddadc8407ffd460e24218e37ef52dd9a3"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.5.16"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.0"
+
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
--- a/test/Project.toml	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/Project.toml	Sat May 20 14:33:25 2023 +0200
@@ -2,6 +2,7 @@
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -10,8 +10,7 @@
 
 @testset "BoundaryOperator" begin
     closure_stencil = Stencil(2.,1.,3.; center = 1)
-    g_1D = EquidistantGrid(11, 0.0, 1.0)
-    g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
+    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
@@ -30,7 +29,7 @@
     end
 
     @testset "Application" begin
-        v = evalOn(g_1D,x->1+x^2)
+        v = eval_on(g_1D,x->1+x^2)
         u = fill(3.124)
         @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3]
         @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2]
@@ -38,7 +37,7 @@
         @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
         @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
 
-        v = evalOn(g_1D, x->1. +x*im)
+        v = eval_on(g_1D, x->1. +x*im)
         @test (op_l*v)[] isa ComplexF64
 
         u = fill(1. +im)
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -9,27 +9,24 @@
 @testset "boundary_restriction" begin
 	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
 	e_closure = parse_stencil(stencil_set["e"]["closure"])
-    g_1D = EquidistantGrid(11, 0.0, 1.0)
-    g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
+    g_1D = equidistant_grid(11, 0.0, 1.0)
+    g_2D = equidistant_grid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
-            @test e_l == boundary_restriction(g_1D,stencil_set,CartesianBoundary{1,Lower}())
+            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
             @test e_l isa LazyTensor{T,0,1} where T
 
-            e_r = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
-            @test e_r == boundary_restriction(g_1D,stencil_set,CartesianBoundary{1,Upper}())
+            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
             @test e_r isa LazyTensor{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D,e_closure,CartesianBoundary{1,Upper}())
-            @test e_w == boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}())
+            e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}())
             @test e_w isa InflatedTensor
             @test e_w isa LazyTensor{T,1,2} where T
         end
@@ -37,8 +34,8 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(e_closure), boundary_identifiers(g_1D))
-            v = evalOn(g_1D,x->1+x^2)
+            e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(stencil_set), boundary_identifiers(g_1D))
+            v = eval_on(g_1D,x->1+x^2)
             u = fill(3.124)
 
             @test (e_l*v)[] == v[1]
@@ -47,7 +44,7 @@
         end
 
         @testset "2D" begin
-            e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(e_closure), boundary_identifiers(g_2D))
+            e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
             v = rand(11, 15)
             u = fill(3.124)
 
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -7,24 +7,23 @@
 import Sbplib.SbpOperators.BoundaryOperator
 
 @testset "normal_derivative" begin
-    g_1D = EquidistantGrid(11, 0.0, 1.0)
-    g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
+    g_1D = equidistant_grid(11, 0.0, 1.0)
+    g_2D = equidistant_grid((11,12), (0.0, 0.0), (1.0,1.0))
     @testset "normal_derivative" begin
     	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    	d_closure = parse_stencil(stencil_set["d1"]["closure"])
         @testset "1D" begin
-            d_l = normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}())
-            @test d_l == normal_derivative(g_1D, stencil_set, CartesianBoundary{1,Lower}())
+            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
             @test d_l isa LazyTensor{T,0,1} where T
         end
         @testset "2D" begin
-            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
-            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
+            d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}())
             Ix = IdentityTensor{Float64}((size(g_2D)[1],))
             Iy = IdentityTensor{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(restrict(g_2D,1),d_closure,CartesianBoundary{1,Lower}())
-            d_r = normal_derivative(restrict(g_2D,2),d_closure,CartesianBoundary{1,Upper}())
+            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}())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
@@ -33,14 +32,13 @@
         end
     end
     @testset "Accuracy" begin
-        v = evalOn(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y)
-        v∂x = evalOn(g_2D, (x,y)-> 2*x + y)
-        v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x)
+        v = eval_on(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y)
+        v∂x = eval_on(g_2D, (x,y)-> 2*x + y)
+        v∂y = eval_on(g_2D, (x,y)-> 2*(y-1) + x)
         # TODO: Test for higher order polynomials?
         @testset "2nd order" begin
         	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-        	d_closure = parse_stencil(stencil_set["d1"]["closure"])
-            d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(d_closure), boundary_identifiers(g_2D))
+            d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
 
             @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ v∂x[end,:] atol = 1e-13
@@ -50,8 +48,7 @@
 
         @testset "4th order" begin
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        	d_closure = parse_stencil(stencil_set["d1"]["closure"])
-            d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(d_closure), boundary_identifiers(g_2D))
+            d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
             
             @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ v∂x[end,:] atol = 1e-13
--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -33,7 +33,7 @@
     @test_throws DomainError ConstantInteriorScalingOperator(4,(2,3), 3)
 
     @testset "Grid constructor" begin
-        g = EquidistantGrid(11, 0., 2.)
+        g = equidistant_grid(11, 0., 2.)
         @test ConstantInteriorScalingOperator(g, 3., (.1,.2)) isa ConstantInteriorScalingOperator{Float64}
     end
 end
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -27,7 +27,7 @@
 end
 
 @testset "undivided_skewed04" begin
-    g = EquidistantGrid(20, 0., 11.)
+    g = equidistant_grid(20, 0., 11.)
     D,Dᵀ = undivided_skewed04(g, 1)
 
     @test D isa LazyTensor{Float64,1,1}
@@ -35,14 +35,14 @@
 
      @testset "Accuracy conditions" begin
         N = 20
-        g = EquidistantGrid(N, 0//1,2//1)
+        g = equidistant_grid(N, 0//1,2//1)
         h = only(spacing(g))
         @testset "D_$p" for p ∈ [1,2,3,4]
             D,Dᵀ = undivided_skewed04(g, p)
 
             @testset "x^$k" for k ∈ 0:p
-                v  = evalOn(g, x->monomial(x,k))
-                vₚₓ = evalOn(g, x->monomial(x,k-p))
+                v  = eval_on(g, x->monomial(x,k))
+                vₚₓ = eval_on(g, x->monomial(x,k-p))
 
                 @test D*v == h^p * vₚₓ
             end
@@ -67,7 +67,7 @@
             return Dmat
         end
 
-        g = EquidistantGrid(11, 0., 1.)
+        g = equidistant_grid(11, 0., 1.)
         @testset "D_$p" for p ∈ [1,2,3,4]
             D,Dᵀ = undivided_skewed04(g, p)
 
@@ -77,6 +77,19 @@
             @test D̄ == D̄ᵀ'
         end
     end
+
+    @testset "2D" begin
+        N = 20
+        g = equidistant_grid((N,2N), (0,0), (2,1))
+        h = spacing.(g.grids)
+
+        D,Dᵀ = undivided_skewed04(g, 3, 2)
+
+        v = eval_on(g, x->monomial(x[1],4)*monomial(x[2],3))
+        d³vdy³ = eval_on(g, x->monomial(x[1],4)*monomial(x[2],0))
+
+        @test D*v ≈ h[2]^3*d³vdy³
+    end
 end
 
 @testset "dissipation_interior_weights" begin
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -24,73 +24,71 @@
     @testset "Constructors" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
 
-        g₁ = EquidistantGrid(11, 0., 1.)
-        g₂ = EquidistantGrid((11,14), (0.,1.), (1.,3.))
+        g₁ = equidistant_grid(11, 0., 1.)
+        g₂ = equidistant_grid((11,14), (0.,1.), (1.,3.))
         
-        @test first_derivative(g₁, stencil_set, 1) isa LazyTensor{Float64,1,1}
+        @test first_derivative(g₁, stencil_set) isa LazyTensor{Float64,1,1}
         @test first_derivative(g₂, stencil_set, 2) isa LazyTensor{Float64,2,2}
-        @test first_derivative(g₁, stencil_set, 1) == first_derivative(g₁, stencil_set)
 
         interior_stencil = CenteredStencil(-1,0,1)
         closure_stencils = [Stencil(-1,1, center=1)]
 
-        @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa LazyTensor{Float64,1,1}
-        @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa VolumeOperator
-        @test first_derivative(g₁, interior_stencil, closure_stencils, 1) == first_derivative(g₁, interior_stencil, closure_stencils)
-        @test first_derivative(g₂, interior_stencil, closure_stencils, 2) isa LazyTensor{Float64,2,2}
+        @test first_derivative(g₁, interior_stencil, closure_stencils) isa LazyTensor{Float64,1,1}
     end
 
     @testset "Accuracy conditions" begin
         N = 20
-        g = EquidistantGrid(N, 0//1,2//1)
+        g = equidistant_grid(N, 0//1,2//1)
         @testset for order ∈ [2,4]
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
-            D₁ = first_derivative(g, stencil_set, 1)
+            D₁ = first_derivative(g, stencil_set)
 
             @testset "boundary x^$k" for k ∈ 0:order÷2
-                v = evalOn(g, x->monomial(x,k))
+                v = eval_on(g, x->monomial(x,k))
 
                 @testset for i ∈ 1:closure_size(D₁)
-                    x, = points(g)[i]
+                    x, = g[i]
                     @test (D₁*v)[i] == monomial(x,k-1)
                 end
 
                 @testset for i ∈ (N-closure_size(D₁)+1):N
-                    x, = points(g)[i]
+                    x, = g[i]
                     @test (D₁*v)[i] == monomial(x,k-1)
                 end
             end
 
             @testset "interior x^$k" for k ∈ 0:order
-                v = evalOn(g, x->monomial(x,k))
+                v = eval_on(g, x->monomial(x,k))
 
-                x, = points(g)[10]
+                x, = g[10]
                 @test (D₁*v)[10] == monomial(x,k-1)
             end
         end
     end
 
     @testset "Accuracy on function" begin
-        # 1D
-        g = EquidistantGrid(30, 0.,1.)
-        v = evalOn(g, x->exp(x))
-        @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)]
-            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
-            D₁ = first_derivative(g, stencil_set, 1)
+        @testset "1D" begin
+            g = equidistant_grid(30, 0.,1.)
+            v = eval_on(g, x->exp(x))
+            @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)]
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
+                D₁ = first_derivative(g, stencil_set)
 
-            @test D₁*v ≈ v rtol=tol
+                @test D₁*v ≈ v rtol=tol
+            end
         end
 
-        # 2D
-        g = EquidistantGrid((30,60), (0.,0.),(1.,2.))
-        v = evalOn(g, (x,y)->exp(0.8x+1.2*y))
-        @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)]
-            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
-            Dx = first_derivative(g, stencil_set, 1)
-            Dy = first_derivative(g, stencil_set, 2)
+        @testset "2D" begin
+            g = equidistant_grid((30,60), (0.,0.),(1.,2.))
+            v = eval_on(g, (x,y)->exp(0.8x+1.2*y))
+            @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)]
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
+                Dx = first_derivative(g, stencil_set, 1)
+                Dy = first_derivative(g, stencil_set, 2)
 
-            @test Dx*v ≈ 0.8v rtol=tol
-            @test Dy*v ≈ 1.2v rtol=tol
+                @test Dx*v ≈ 0.8v rtol=tol
+                @test Dy*v ≈ 1.2v rtol=tol
+            end
         end
     end
 end
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -15,24 +15,18 @@
     closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
     Lx = 3.5
     Ly = 3.
-    g_1D = EquidistantGrid(121, 0.0, Lx)
-    g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
+    g_1D = equidistant_grid(121, 0.0, Lx)
+    g_2D = equidistant_grid((121,123), (0.0, 0.0), (Lx, Ly))
 
     @testset "Constructors" begin
         @testset "1D" begin
-            Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils,1)
-            @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils)
-            @test Dₓₓ == second_derivative(g_1D,stencil_set,1)
-            @test Dₓₓ == second_derivative(g_1D,stencil_set)
-            @test Dₓₓ isa VolumeOperator
+            Dₓₓ = second_derivative(g_1D, stencil_set)
+            @test Dₓₓ == second_derivative(g_1D, inner_stencil, closure_stencils)
+            @test Dₓₓ isa LazyTensor{Float64,1,1}
         end
         @testset "2D" begin
-            Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1)
-            D2 = second_derivative(g_1D,inner_stencil,closure_stencils,1)
-            I = IdentityTensor{Float64}(size(g_2D)[2])
-            @test Dₓₓ == D2⊗I
-            @test Dₓₓ == second_derivative(g_2D,stencil_set,1)
-            @test Dₓₓ isa LazyTensor{T,2,2} where T
+            Dₓₓ = second_derivative(g_2D,stencil_set,1)
+            @test Dₓₓ isa LazyTensor{Float64,2,2}
         end
     end
 
@@ -45,10 +39,10 @@
             maxOrder = 4;
             for i = 0:maxOrder-1
                 f_i(x) = 1/factorial(i)*x^i
-                monomials = (monomials...,evalOn(g_1D,f_i))
+                monomials = (monomials...,eval_on(g_1D,f_i))
             end
-            v = evalOn(g_1D,x -> sin(x))
-            vₓₓ = evalOn(g_1D,x -> -sin(x))
+            v = eval_on(g_1D,x -> sin(x))
+            vₓₓ = eval_on(g_1D,x -> -sin(x))
 
             # 2nd order interior stencil, 1nd order boundary stencil,
             # implies that L*v should be exact for monomials up to order 2.
@@ -77,15 +71,15 @@
         end
 
         @testset "2D" begin
-            l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
+            l2(v) = sqrt(prod(spacing.(g_2D.grids))*sum(v.^2));
             binomials = ()
             maxOrder = 4;
             for i = 0:maxOrder-1
                 f_i(x,y) = 1/factorial(i)*y^i + x^i
-                binomials = (binomials...,evalOn(g_2D,f_i))
+                binomials = (binomials...,eval_on(g_2D,f_i))
             end
-            v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
-            v_yy = evalOn(g_2D,(x,y) -> -cos(y))
+            v = eval_on(g_2D, (x,y) -> sin(x)+cos(y))
+            v_yy = eval_on(g_2D,(x,y) -> -cos(y))
 
             # 2nd order interior stencil, 1st order boundary stencil,
             # implies that L*v should be exact for binomials up to order 2.
@@ -94,7 +88,7 @@
                 Dyy = second_derivative(g_2D,stencil_set,2)
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
                 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
-                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
+                @test Dyy*binomials[3] ≈ eval_on(g_2D,(x,y)->1.) atol = 5e-9
                 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
             end
 
@@ -107,8 +101,8 @@
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
                 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
-                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
-                @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
+                @test Dyy*binomials[3] ≈ eval_on(g_2D,(x,y)->1.) atol = 5e-9
+                @test Dyy*binomials[4] ≈ eval_on(g_2D,(x,y)->y) atol = 5e-9
                 @test Dyy*v ≈ v_yy rtol = 5e-4 norm = l2
             end
         end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -10,47 +10,39 @@
     Lx = π/2.
     Ly = Float64(π)
     Lz = 1.
-    g_1D = EquidistantGrid(77, 0.0, Lx)
-    g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-    g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
-    integral(H,v) = sum(H*v)
+    g_1D = equidistant_grid(77, 0.0, Lx)
+    g_2D = equidistant_grid((77,66), (0.0, 0.0), (Lx,Ly))
+    g_3D = equidistant_grid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
     @testset "inner_product" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-        quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
         @testset "0D" begin
-            H = inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
-            @test H == inner_product(EquidistantGrid{Float64}(), stencil_set)
-            @test H == IdentityTensor{Float64}()
+            H = inner_product(ZeroDimGrid(0.), stencil_set)
             @test H isa LazyTensor{T,0,0} where T
         end
         @testset "1D" begin
-            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
-            @test H == inner_product(g_1D, stencil_set)
-            @test H isa ConstantInteriorScalingOperator
+            H = inner_product(g_1D, stencil_set)
             @test H isa LazyTensor{T,1,1} where T
         end
         @testset "2D" begin
-            H = inner_product(g_2D, quadrature_interior, quadrature_closure)
-            H_x = inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure)
-            H_y = inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure)
-            @test H == inner_product(g_2D, stencil_set)
+            H = inner_product(g_2D, stencil_set)
+            H_x = inner_product(g_2D.grids[1], stencil_set)
+            H_y = inner_product(g_2D.grids[2], stencil_set)
             @test H == H_x⊗H_y
             @test H isa LazyTensor{T,2,2} where T
         end
+
+        # TBD: Should there be more tests?
     end
 
     @testset "Sizes" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-        quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
         @testset "1D" begin
-            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+            H = inner_product(g_1D, stencil_set)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = inner_product(g_2D, quadrature_interior, quadrature_closure)
+            H = inner_product(g_2D, stencil_set)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -61,52 +53,44 @@
             v = ()
             for i = 0:4
                 f_i(x) = 1/factorial(i)*x^i
-                v = (v...,evalOn(g_1D,f_i))
+                v = (v...,eval_on(g_1D,f_i))
             end
-            u = evalOn(g_1D,x->sin(x))
+            u = eval_on(g_1D,x->sin(x))
 
             @testset "2nd order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+                H = inner_product(g_1D, stencil_set)
                 for i = 1:2
-                    @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
+                    @test sum(H*v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
-                @test integral(H,u) ≈ 1. rtol = 1e-4
+                @test sum(H*u) ≈ 1. rtol = 1e-4
             end
 
             @testset "4th order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+                H = inner_product(g_1D, stencil_set)
                 for i = 1:4
-                    @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
+                    @test sum(H*v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
-                @test integral(H,u) ≈ 1. rtol = 1e-8
+                @test sum(H*u) ≈ 1. rtol = 1e-8
             end
         end
 
         @testset "2D" begin
             b = 2.1
             v = b*ones(Float64, size(g_2D))
-            u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
+            u = eval_on(g_2D,(x,y)->sin(x)+cos(y))
             @testset "2nd order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
-                @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
-                @test integral(H,u) ≈ π rtol = 1e-4
+                H = inner_product(g_2D, stencil_set)
+                @test sum(H*v) ≈ b*Lx*Ly rtol = 1e-13
+                @test sum(H*u) ≈ π rtol = 1e-4
             end
             @testset "4th order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
-                @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
-                @test integral(H,u) ≈ π rtol = 1e-8
+                H = inner_product(g_2D, stencil_set)
+                @test sum(H*v) ≈ b*Lx*Ly rtol = 1e-13
+                @test sum(H*u) ≈ π rtol = 1e-8
             end
         end
     end
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -9,29 +9,22 @@
 @testset "Diagonal-stencil inverse_inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
-    g_1D = EquidistantGrid(77, 0.0, Lx)
-    g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    g_1D = equidistant_grid(77, 0.0, Lx)
+    g_2D = equidistant_grid((77,66), (0.0, 0.0), (Lx,Ly))
     @testset "inverse_inner_product" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-        quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
         @testset "0D" begin
-            Hi = inverse_inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
-            @test Hi == inverse_inner_product(EquidistantGrid{Float64}(), stencil_set)
-            @test Hi == IdentityTensor{Float64}()
+            Hi = inverse_inner_product(ZeroDimGrid(1.), stencil_set)
             @test Hi isa LazyTensor{T,0,0} where T
         end
         @testset "1D" begin
-            Hi = inverse_inner_product(g_1D,  quadrature_interior, quadrature_closure)
-            @test Hi == inverse_inner_product(g_1D, stencil_set)
-            @test Hi isa ConstantInteriorScalingOperator
+            Hi = inverse_inner_product(g_1D, stencil_set)
             @test Hi isa LazyTensor{T,1,1} where T
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
-            Hi_x = inverse_inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure)
-            Hi_y = inverse_inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure)
-            @test Hi == inverse_inner_product(g_2D, stencil_set)
+            Hi = inverse_inner_product(g_2D, stencil_set)
+            Hi_x = inverse_inner_product(g_2D.grids[1], stencil_set)
+            Hi_y = inverse_inner_product(g_2D.grids[2], stencil_set)
             @test Hi == Hi_x⊗Hi_y
             @test Hi isa LazyTensor{T,2,2} where T
         end
@@ -39,15 +32,13 @@
 
     @testset "Sizes" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-        quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
         @testset "1D" begin
-            Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
+            Hi = inverse_inner_product(g_1D, stencil_set)
             @test domain_size(Hi) == size(g_1D)
             @test range_size(Hi) == size(g_1D)
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
+            Hi = inverse_inner_product(g_2D, stencil_set)
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -55,45 +46,37 @@
 
     @testset "Accuracy" begin
         @testset "1D" begin
-            v = evalOn(g_1D,x->sin(x))
-            u = evalOn(g_1D,x->x^3-x^2+1)
+            v = eval_on(g_1D,x->sin(x))
+            u = eval_on(g_1D,x->x^3-x^2+1)
             @testset "2nd order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
-                Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
+                H = inner_product(g_1D, stencil_set)
+                Hi = inverse_inner_product(g_1D, stencil_set)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
-                Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
+                H = inner_product(g_1D, stencil_set)
+                Hi = inverse_inner_product(g_1D, stencil_set)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
         end
         @testset "2D" begin
-            v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
-            u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
+            v = eval_on(g_2D,(x,y)->sin(x)+cos(y))
+            u = eval_on(g_2D,(x,y)->x*y + x^5 - sqrt(y))
             @testset "2nd order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
-                Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
+                H = inner_product(g_2D, stencil_set)
+                Hi = inverse_inner_product(g_2D, stencil_set)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
-                quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
-                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
-                Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
+                H = inner_product(g_2D, stencil_set)
+                Hi = inverse_inner_product(g_2D, stencil_set)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -4,40 +4,40 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
-# Default stencils (4th order)
-operator_path = sbp_operators_path()*"standard_diagonal.toml"
-stencil_set = read_stencil_set(operator_path; order=4)
-inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-g_1D = EquidistantGrid(101, 0.0, 1.)
-g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+@testset "Laplace" begin
+    # Default stencils (4th order)
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    stencil_set = read_stencil_set(operator_path; order=4)
+    g_1D = equidistant_grid(101, 0.0, 1.)
+    g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
 
-@testset "Laplace" begin
     @testset "Constructors" begin
         @testset "1D" begin
-            Δ = laplace(g_1D, inner_stencil, closure_stencils)
-            @test Laplace(g_1D, stencil_set) == Laplace(Δ, stencil_set)
-            @test Laplace(g_1D, stencil_set) isa LazyTensor{T,1,1}  where T
+            @test Laplace(g_1D, stencil_set) == Laplace(laplace(g_1D, stencil_set), stencil_set)
+            @test Laplace(g_1D, stencil_set) isa LazyTensor{Float64,1,1}
         end
         @testset "3D" begin
-            Δ = laplace(g_3D, inner_stencil, closure_stencils)
-            @test Laplace(g_3D, stencil_set) == Laplace(Δ,stencil_set)
-            @test Laplace(g_3D, stencil_set) isa LazyTensor{T,3,3} where T
+            @test Laplace(g_3D, stencil_set) == Laplace(laplace(g_3D, stencil_set),stencil_set)
+            @test Laplace(g_3D, stencil_set) isa LazyTensor{Float64,3,3}
         end
     end
 
     # Exact differentiation is measured point-wise. In other cases
     # the error is measured in the l2-norm.
     @testset "Accuracy" begin
-        l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2));
+        l2(v) = sqrt(prod(spacing.(g_3D.grids))*sum(v.^2));
         polynomials = ()
         maxOrder = 4;
         for i = 0:maxOrder-1
             f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i)
-            polynomials = (polynomials...,evalOn(g_3D,f_i))
+            polynomials = (polynomials...,eval_on(g_3D,f_i))
         end
-        v = evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
-        Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
+        # v = eval_on(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
+        # Δv = eval_on(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
+
+        v =  eval_on(g_3D, x̄ -> sin(x̄[1]) + cos(x̄[2]) + exp(x̄[3]))
+        Δv = eval_on(g_3D, x̄ -> -sin(x̄[1]) - cos(x̄[2]) + exp(x̄[3]))
+        @inferred v[1,2,3]
 
         # 2nd order interior stencil, 1st order boundary stencil,
         # implies that L*v should be exact for binomials up to order 2.
@@ -67,19 +67,24 @@
 end
 
 @testset "laplace" begin
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    stencil_set = read_stencil_set(operator_path; order=4)
+    g_1D = equidistant_grid(101, 0.0, 1.)
+    g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+
     @testset "1D" begin
-        Δ = laplace(g_1D, inner_stencil, closure_stencils)
-        @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils, 1)
-        @test Δ isa LazyTensor{T,1,1}  where T
+        Δ = laplace(g_1D, stencil_set)
+        @test Δ == second_derivative(g_1D, stencil_set)
+        @test Δ isa LazyTensor{Float64,1,1}
     end
     @testset "3D" begin
-        Δ = laplace(g_3D, inner_stencil, closure_stencils)
-        @test Δ isa LazyTensor{T,3,3} where T
-        Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
-        Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
-        Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
+        Δ = laplace(g_3D, stencil_set)
+        @test Δ isa LazyTensor{Float64,3,3}
+        Dxx = second_derivative(g_3D, stencil_set, 1)
+        Dyy = second_derivative(g_3D, stencil_set, 2)
+        Dzz = second_derivative(g_3D, stencil_set, 3)
         @test Δ == Dxx + Dyy + Dzz
-        @test Δ isa LazyTensor{T,3,3} where T
+        @test Δ isa LazyTensor{Float64,3,3}
     end
 end
 
--- a/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -6,40 +6,9 @@
 
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.StencilOperatorDistinctClosures
-import Sbplib.SbpOperators.stencil_operator_distinct_closures
-
-@testset "stencil_operator_distinct_closures" begin
-    lower_closure = (
-        Stencil(-1,1, center=1),
-    )
-
-    inner_stencil = Stencil(-2,2, center=1)
-
-    upper_closure = (
-        Stencil(-3,3, center=1),
-        Stencil(-4,4, center=2),
-    )
-
-    g₁ = EquidistantGrid(5, 0., 1.)
-    g₂ = EquidistantGrid((5,5), (0.,0.), (1.,1.))
-    h = 1/4
-
-    A₁  = stencil_operator_distinct_closures(g₁, inner_stencil, lower_closure, upper_closure, 1)
-    A₂¹ = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 1)
-    A₂² = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 2)
-
-    v₁ = evalOn(g₁, x->x)
-
-    u = [1., 2., 2., 3., 4.]*h
-    @test A₁*v₁ == u
-
-    v₂ = evalOn(g₂, (x,y)-> x + 3y)
-    @test A₂¹*v₂ == repeat(u, 1, 5)
-    @test A₂²*v₂ == repeat(3u', 5, 1)
-end
 
 @testset "StencilOperatorDistinctClosures" begin
-    g = EquidistantGrid(11, 0., 1.)
+    g = equidistant_grid(11, 0., 1.)
 
     lower_closure = (
         Stencil(-1,1, center=1),
--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Sat May 20 14:33:25 2023 +0200
@@ -14,7 +14,7 @@
 @testset "VolumeOperator" begin
     inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
     closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2))
-    g = EquidistantGrid(11,0.,1.)
+    g = equidistant_grid(11,0.,1.)
 
     @testset "Constructors" begin
         op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)
--- a/test/runtests.jl	Tue Apr 04 21:46:06 2023 +0200
+++ b/test/runtests.jl	Sat May 20 14:33:25 2023 +0200
@@ -49,4 +49,5 @@
 
 @testset "$testsetname" begin
     run_testfiles(ARGS)
+    println()
 end