--- a/Manifest.toml	Mon Nov 30 16:28:32 2020 +0100
+++ b/Manifest.toml	Fri Jan 01 16:34:55 2021 +0100
@@ -1,12 +1,29 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
 [[OffsetArrays]]
 git-tree-sha1 = "9011c7c98769c451f83869a4d66461e2f23bc80b"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 version = "1.2.1"
 
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[TOML]]
+deps = ["Dates"]
+git-tree-sha1 = "d0ac7eaad0fb9f6ba023a1d743edca974ae637c4"
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
+
 [[TiledIteration]]
 deps = ["OffsetArrays"]
 git-tree-sha1 = "98693daea9bb49aba71eaad6b168b152d2310358"
 uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 version = "0.2.4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

--- a/Notes.md	Mon Nov 30 16:28:32 2020 +0100
+++ b/Notes.md	Fri Jan 01 16:34:55 2021 +0100
@@ -11,23 +11,31 @@
  Maybe if we should have dynamic sizing it could be only for the range. `domain_size` would not be implemented. And the `range_size` would be a function of a vector that the TensorMapping is applied to.
 
 ## Reasearch and thinking
- - [ ] Use a trait to indicate if a TensorMapping uses indices with regions.
-    The default should be that they do NOT.
-        - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle?
  - [ ] Use a trait to indicate that a TensorMapping har the same range and domain?
  - [ ] Rename all the Tensor stuff to just LazyOperator, LazyApplication and so on?
- - [ ] Figure out repeated application of regioned TensorMappings. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work?
  - [ ] Check how the native julia doc generator works
     - [ ] Check if Vidars design docs fit in there
  - [ ] Create a macro @lazy which replaces a binary op (+,-) by its lazy equivalent? Would be a neat way to indicate which evaluations are lazy without cluttering/confusing with special characters.
  - [ ] Specificera operatorer i TOML eller något liknande?
  H.. H_gamma etc.)
- - [ ] Dispatch in Lower() instead of the type Lower so `::Lower` instead of `::Type{Lower}` ???
+ - [ ] 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 TensorMapping 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?
 
+## Regions and tensormappings
+- [ ] Use a trait to indicate if a TensorMapping uses indices with regions.
+    The default should be that they do NOT.
+        - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle?
+ - [ ] Figure out repeated application of regioned TensorMappings. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work?
+
+## Boundschecking and dimension checking
+Does it make sense to have boundschecking only in getindex methods?
+This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway.
+
+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:
 
@@ -101,34 +109,50 @@
 ### Test-applikationer
 div och grad operationer
 
-### Alternativ:
+Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensor fä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?
 
-#### 1.Använda tuplar
-Fördelar:
+Ä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?
 
-Nackdelar:
+### 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.
+
+En nackdel kan vara hur man ska få ut gridfunktionen för tex andra komponenten.
 
 Syntax:
 ```
-f(x,y) = x^2 + y^2
-gf = eval_on_grid(g,f)
-gf[2,3] # En tupel för en given gridpunkt
-gf[2,3][2] # Andra komponenten av vektor fältet i en punkt.
-
+f(x̄) = x̄
+gf = evalOn(g, f)
+gf[2,3] # x̄ för en viss gridpunkt
+gf[2,3][2] # x̄[2] för en viss gridpunkt
 ```
 
-#### 2.AbstractArray{T,2} where T
-Låta alla saker ta in AbstractArray{T,2} where T. Där T kan vara lite vad som helst, tillexemel en SVector eller Array. Men Differens-opertorerna bryr sig inte om det.
+Note: Behöver bestämma om eval on skickar in `x̄` eller `x̄...` till `f`. Eller om man kan stödja båda.
 
-En nackdel kan var hur man ska få ut gridfunktionen för tex andra komponenten.
+### 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.
+
+TBD: Just nu gör `apply_transpose` antagandet att domän-typen är samma som range-typen. Det behöver vi på något sätt bryta. Ett alternativ är låta en TensorMapping ha `T_domain` och `T_range` istället för bara `T`. Känns dock lite grötigt. Ett annat alternativ skulle vara någon typ av trait för transpose? Den skulle kunna innehålla typen som transponatet agerar på? Vet inte om det fungerar dock.
+
+TBD: Vad är målet med `T`-parametern för en TensorMapping? Om vi vill kunna applicera en difference operator på vad som helst kan man inte anta att en `TensorMapping{T}` bara agerar på instanser av `T`.
 
-Syntax:
-```
-gf = eval(...)
-gf[2,3] # Hela vektorn för en gridpunkt
-gf[2,3][2] # Andra komponenten av vektor fältet.
-```
-#### 3.AbstractArray{T,2+1} where T
+Man kan implementera `∇` som en tensormapping som agerar på T och returnerar `StaticVector{N,T} where N`.
+(Man skulle eventuellt också kunna låta den agera på `StaticMatrix{N,T,D} where N` och returnera `StaticMatrix{M,T,D+1}`. Frågan är om man vinner något på det...)
+
+Skulle kunna ha en funktion `range_type(::TensorMapping, ::Type{domain_type})`
+
+Kanske kan man implementera `⋅(tm::TensorMapping{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:
@@ -150,7 +174,7 @@
 Påverkas detta av hur vi förväntar oss kunna skapa lata gridfunktioner?
 
 ### Komponenter som gridfunktioner
-För alternativ 1 och 2 har vi problemet hur vi får ut komponenter som vektorfält. Detta behöver antagligen kunna ske lazy.
+En viktig operation för vektor fält är att kunna få ut komponenter som grid-funktioner. Detta behöver antagligen kunna ske lazy.
 Det finns ett par olika lösningar:
 * Implementera en egen typ av view som tar hand om detta. Eller Accessors.jl?
 * Använda en TensorMapping

--- a/Project.toml	Mon Nov 30 16:28:32 2020 +0100
+++ b/Project.toml	Fri Jan 01 16:34:55 2021 +0100
@@ -4,6 +4,7 @@
 version = "0.1.0"
 
 [deps]
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
 [compat]

--- a/TODO.md	Mon Nov 30 16:28:32 2020 +0100
+++ b/TODO.md	Fri Jan 01 16:34:55 2021 +0100
@@ -9,11 +9,18 @@
  - [ ] Add 1D operators (D1, D2, e, d ... ) as TensorOperators
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
  - [ ] Add a quick and simple way of running all tests for all subpackages.
- - [ ] Replace getindex hack for flatteing tuples with flatten_tuple.
+ - [ ] Replace getindex hack for flatteing tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`)
  - [ ] Use `@inferred` in a lot of tests.
  - [ ] Make sure we are setting tolerances in tests in a consistent way
  - [ ] Add check for correct domain sizes to lazy tensor operations using SizeMismatch
- - [ ] Write down some coding guideline or checklist for code convetions. For example i,j,... för indecies and I for multi-index
+ - [ ] Write down some coding guideline or checklist for code convetions. For example i,j,... for indecies and I for multi-index
+ - [ ] Add boundschecking in TensorMappingApplication
+ - [ ] Start renaming things in LazyTensors
+ - [ ] Clean up RegionIndices
+    1. [ ] Write tests for how things should work
+    2. [ ] Update RegionIndices accordingly
+    3. [ ] Fix the rest of the library
+ - [ ] Add posibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`?
 
 ## Repo
  - [ ] Add Vidar to the authors list

--- a/src/Grids/Grids.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/src/Grids/Grids.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -7,7 +7,7 @@
 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  #TODO: Should return R()
+region(::CartesianBoundary{Dim, R}) where {Dim, R} = R()
 
 export dim, region
 

--- a/src/SbpOperators/BoundaryValue.jl	Mon Nov 30 16:28:32 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,127 +0,0 @@
-"""
-    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
-
-Implements the boundary operator `e` as a TensorMapping
-"""
-struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
-    eClosure::Stencil{T,M}
-    bId::CartesianBoundary
-end
-export BoundaryValue
-
-# TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
-# Can we give special treatment to TensorMappings that go to a higher dim?
-function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T
-    if dim(e.bId) == 1
-        return (UnknownDim, domain_size[1])
-    elseif dim(e.bId) == 2
-        return (domain_size[1], UnknownDim)
-    end
-end
-LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
-# TODO: Make a nicer solution for 3-dim(e.bId)
-
-# TODO: Make this independent of dimension
-function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
-    i = I[dim(e.bId)]
-    j = I[3-dim(e.bId)]
-    N_i = size(e.grid)[dim(e.bId)]
-    return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
-end
-
-function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
-    u = selectdim(v,3-dim(e.bId),Int(I[1]))
-    return apply_boundary_value_transpose(e.op, u, region(e.bId))
-end
-
-function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    apply_stencil(op.eClosure,v,1)
-end
-
-function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    apply_stencil_backwards(op.eClosure,v,length(v))
-end
-export apply_boundary_value_transpose
-
-function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(Int(i)) <= N)
-        throw(BoundsError())
-    end
-    op.eClosure[Int(i)-1]*v
-end
-
-function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index,  N::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(Int(i)) <= N)
-        throw(BoundsError())
-    end
-    op.eClosure[N-Int(i)]*v
-end
-export apply_boundary_value
-
-
-"""
-    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
-
-Implements the boundary operator `e` as a TensorMapping
-"""
-struct BoundaryValue{D,T,M,R} <: TensorMapping{T,D,1}
-    e:BoundaryOperator{T,M,R}
-    bId::CartesianBoundary
-end
-
-function LazyTensors.apply_transpose(bv::BoundaryValue{T,M,Lower}, v::AbstractVector{T}, i::Index) where T
-    u = selectdim(v,3-dim(bv.bId),Int(I[1]))
-    return apply_transpose(bv.e, u, I)
-end
-
-
-"""
-    BoundaryOperator{T,N,R} <: TensorMapping{T,1,1}
-
-Implements the boundary operator `e` as a TensorMapping
-"""
-export BoundaryOperator
-struct BoundaryOperator{T,M,R<:Region} <: TensorMapping{T,1,1}
-    closure::Stencil{T,M}
-end
-
-function LazyTensors.range_size(e::BoundaryOperator, domain_size::NTuple{1,Integer})
-    return UnknownDim
-end
-
-LazyTensors.domain_size(e::BoundaryOperator{T}, range_size::NTuple{1,Integer}) where T = range_size
-
-function LazyTensors.apply_transpose(e::BoundaryOperator{T,M,Lower}, v::AbstractVector{T}, i::Index{Lower}) where T
-    @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here?
-        throw(BoundsError())
-    end
-    apply_stencil(e.closure,v,Int(i))
-end
-
-function LazyTensors.apply_transpose(e::BoundaryOperator{T,M,Upper}}, v::AbstractVector{T}, i::Index{Upper}) where T
-    @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here?
-        throw(BoundsError())
-    end
-    apply_stencil_backwards(e.closure,v,Int(i))
-end
-
-function LazyTensors.apply_transpose(e::BoundaryOperator{T}, v::AbstractVector{T}, i::Index) where T
-    @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here?
-        throw(BoundsError())
-    end
-    return eltype(v)(0)
-end
-
-#TODO: Implement apply in a meaningful way. Should it return a vector or a single value (perferable?) Should fit into  the
-function LazyTensors.apply(e::BoundaryOperator, v::AbstractVector, i::Index)
-    @boundscheck if !(0<length(Int(i)) <= length(v))
-        throw(BoundsError())
-    end
-    return e.closure[Int(i)].*v
-end

--- a/src/SbpOperators/SbpOperators.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -12,5 +12,6 @@
 include("laplace/laplace.jl")
 include("quadrature/diagonal_quadrature.jl")
 include("quadrature/inverse_diagonal_quadrature.jl")
+include("boundaryops/boundary_restriction.jl")
 
 end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -0,0 +1,81 @@
+"""
+    boundary_restriction(grid,closureStencil,boundary)
+
+Creates a boundary restriction operator on a `Dim`-dimensional grid for the
+specified `boundary`.
+
+When `Dim=1`, the corresponding `BoundaryRestriction` tensor mapping is returned.
+When `Dim>1`, the `BoundaryRestriction` `e` is inflated by the outer product
+of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
+the boundary restriction operator in the y-direction direction is `Ix⊗e⊗Iz`.
+"""
+function boundary_restriction(grid::EquidistantGrid{Dim,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary) where {Dim,T,M}
+    # Create 1D boundary restriction operator
+    r = region(boundary)
+    d = dim(boundary)
+    e = BoundaryRestriction(restrict(grid, d), closureStencil, r)
+
+    # Create 1D IdentityMappings for each coordinate direction
+    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
+    Is = IdentityMapping{T}.(size.(one_d_grids))
+
+    # Formulate the correct outer product sequence of the identity mappings and
+    # the boundary restriction operator
+    parts = Base.setindex(Is, e, d)
+    return foldl(⊗, parts)
+end
+
+export boundary_restriction
+
+"""
+    BoundaryRestriction{T,R,N} <: TensorMapping{T,0,1}
+
+Implements the boundary operator `e` for 1D as a `TensorMapping`
+
+`e` is the restriction of a grid function to the boundary using some `closureStencil`.
+The boundary to restrict to is determined by `R`.
+
+`e'` is the prolongation of a zero dimensional array to the whole grid using the same `closureStencil`.
+"""
+struct BoundaryRestriction{T,R<:Region,N} <: TensorMapping{T,0,1}
+    stencil::Stencil{T,N}
+    size::Int
+end
+export BoundaryRestriction
+
+BoundaryRestriction{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryRestriction{T,R,N}(stencil, size)
+
+function BoundaryRestriction(grid::EquidistantGrid{1}, closureStencil::Stencil{T,N}, region::Region) where {T,N}
+    return BoundaryRestriction{T,typeof(region),N}(closureStencil,size(grid)[1])
+end
+
+closure_size(::BoundaryRestriction{T,R,N}) where {T,R,N} = N
+
+LazyTensors.range_size(e::BoundaryRestriction) = ()
+LazyTensors.domain_size(e::BoundaryRestriction) = (e.size,)
+
+function LazyTensors.apply(e::BoundaryRestriction{T,Lower}, v::AbstractVector{T}) where T
+    apply_stencil(e.stencil,v,1)
+end
+
+function LazyTensors.apply(e::BoundaryRestriction{T,Upper}, v::AbstractVector{T}) where T
+    apply_stencil_backwards(e.stencil,v,e.size)
+end
+
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T
+    return e.stencil[Int(i)-1]*v[]
+end
+
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T
+    return e.stencil[e.size[1] - Int(i)]*v[]
+end
+
+# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods.
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i::Index) where T
+    return zero(T)
+end
+
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i) where T
+    r = getregion(i, closure_size(e), e.size)
+    apply_transpose(e, v, Index(i,r))
+end

--- a/src/SbpOperators/d2.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/src/SbpOperators/d2.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -1,4 +1,4 @@
-export D2, closuresize, readOperator
+export D2, closuresize
 
 @enum Parity begin
     odd = -1

--- a/src/SbpOperators/operators/d2_2nd.txt	Mon Nov 30 16:28:32 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,13 +0,0 @@
-# D2 order 2
-
-boundary_stencils
-1 -2 1
-
-inner_stencil
-1 -2 1
-
-e
-1
-
-d1
--3/2 2 -1/2

--- a/src/SbpOperators/operators/d2_4th.txt	Mon Nov 30 16:28:32 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,16 +0,0 @@
-# D2 order 4
-
-boundary_stencils
-     2    -5       4      -1     0     0
-     1    -2       1       0     0     0
- -4/43 59/43 -110/43   59/43 -4/43     0
- -1/49     0   59/49 -118/49 64/49 -4/49
-
-inner_stencil
--1/12     4/3  -5/2   4/3 -1/12
-
-e
-1
-
-d1
--11/6 3 -3/2 1/3
\ No newline at end of file

--- a/src/SbpOperators/operators/h_2nd.txt	Mon Nov 30 16:28:32 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-# H order 2
-
-closure
-1/2
-
-inner_stencil
-1

--- a/src/SbpOperators/operators/h_4th.txt	Mon Nov 30 16:28:32 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-# H order 4
-
-closure
-17/48 59/48 43/48 49/48
-
-inner_stencil
-1
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Fri Jan 01 16:34:55 2021 +0100
@@ -0,0 +1,36 @@
+[meta]
+authors = "Ken Mattson"
+descripion = "Standard operators for equidistant grids"
+type = "equidistant"
+
+[order2]
+H.inner = ["1"]
+H.closure = ["1/2"]
+
+D1.inner_stencil = ["-1/2", "0", "1/2"]
+D1.closure_stencils = [
+    ["-1", "1"],
+]
+
+D2.inner_stencil = ["1", "-2", "1"]
+D2.closure_stencils = [
+    ["1", "-2", "1"],
+]
+
+e.closure = ["1"]
+d1.closure = ["-3/2", "2", "-1/2"]
+
+[order4]
+H.inner = ["1"]
+H.closure = ["17/48", "59/48", "43/48", "49/48"]
+
+D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
+D2.closure_stencils = [
+    [     "2",    "-5",      "4",       "-1",     "0",     "0"],
+    [     "1",    "-2",      "1",        "0",     "0",     "0"],
+    [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"],
+    [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"],
+]
+
+e.closure = ["1"]
+d1.closure = ["-11/6", "3", "-3/2", "1/3"]

--- a/src/SbpOperators/readoperator.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/src/SbpOperators/readoperator.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -1,30 +1,39 @@
-function readOperator(D2fn, Hfn)
-    d = readSectionedFile(D2fn)
-    h = readSectionedFile(Hfn)
+using TOML
+
+export read_D2_operator
+export read_stencil
+export read_stencils
+export read_tuple
+
+export get_stencil
+export get_stencils
+export get_tuple
+
+
+function read_D2_operator(fn; order)
+    operators = TOML.parsefile(fn)["order$order"]
+    D2 = operators["D2"]
+    H = operators["H"]
+    e = operators["e"]
+    d1 = operators["d1"]
 
     # Create inner stencil
-    innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
-    width = length(innerStencilWeights)
-    r = (-div(width,2), div(width,2))
-
-    innerStencil = Stencil(r, innerStencilWeights)
+    innerStencil = get_stencil(operators, "D2", "inner_stencil")
 
     # Create boundary stencils
-    boundarySize = length(d["boundary_stencils"])
+    boundarySize = length(D2["closure_stencils"])
     closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
 
     for i ∈ 1:boundarySize
-        stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
-        width = length(stencilWeights)
-        r = (1-i,width-i)
-        closureStencils = (closureStencils..., Stencil(r, stencilWeights))
+        closureStencils = (closureStencils..., get_stencil(operators, "D2", "closure_stencils", i; center=i))
     end
 
-    quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
-    eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
-    dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
+    # TODO: Get rid of the padding here. Any padding should be handled by the consturctor accepting the stencils.
+    quadratureClosure = pad_tuple(toml_string_array_to_tuple(Float64, H["closure"]), boundarySize)
+    eClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, e["closure"]), boundarySize), center=1)
+    dClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, d1["closure"]), boundarySize), center=1)
 
-    d2 = D2(
+    d2 = SbpOperators.D2(
         quadratureClosure,
         innerStencil,
         closureStencils,
@@ -36,36 +45,100 @@
     return d2
 end
 
-function readSectionedFile(filename)::Dict{String, Vector{String}}
-    f = open(filename)
-    sections = Dict{String, Vector{String}}()
-    currentKey = ""
+
+"""
+    read_stencil(fn, path...; [center])
+
+Read a stencil at `path` from the file with name `fn`.
+If a center is specified the given element of the stecil is set as the center.
+
+See also: [`read_stencils`](@ref), [`read_tuple`](@ref), [`get_stencil`](@ref).
 
-    for ln ∈ eachline(f)
-        if ln == "" || ln[1] == '#' # Skip comments and empty lines
-            continue
-        end
+# Examples
+```
+read_stencil(sbp_operators_path()*"standard_diagonal.toml", "order2", "D2", "inner_stencil")
+read_stencil(sbp_operators_path()*"standard_diagonal.toml", "order2", "d1", "closure"; center=1)
+```
+"""
+read_stencil(fn, path...; center=nothing) = get_stencil(TOML.parsefile(fn), path...; center=center)
+
+"""
+    read_stencils(fn, path...; centers)
+
+Read stencils at `path` from the file `fn`.
+Centers of the stencils are specified as a tuple or array in `centers`.
 
-        if isletter(ln[1]) # Found start of new section
-            if ~haskey(sections, ln)
-                sections[ln] =  Vector{String}()
-            end
-            currentKey = ln
-            continue
-        end
+See also: [`read_stencil`](@ref), [`read_tuple`](@ref), [`get_stencils`](@ref).
+"""
+read_stencils(fn, path...; centers) = get_stencils(TOML.parsefile(fn), path...; centers=centers)
+
+"""
+    read_tuple(fn, path...)
+
+Read tuple at `path` from the file `fn`.
+
+See also: [`read_stencil`](@ref), [`read_stencils`](@ref), [`get_tuple`](@ref).
+"""
+read_tuple(fn, path...) = get_tuple(TOML.parsefile(fn), path...)
 
-        push!(sections[currentKey], ln)
+"""
+    get_stencil(parsed_toml, path...; center=nothing)
+
+Same as [`read_stencil`](@ref)) but takes already parsed toml.
+"""
+get_stencil(parsed_toml, path...; center=nothing) = get_stencil(parsed_toml[path[1]], path[2:end]...; center=center)
+function get_stencil(parsed_toml; center=nothing)
+    @assert parsed_toml isa Vector{String}
+    stencil_weights = Float64.(parse_rational.(parsed_toml))
+
+    width = length(stencil_weights)
+
+    if isnothing(center)
+        center = div(width,2)+1
     end
 
-    return sections
+    return Stencil(Tuple(stencil_weights), center=center)
 end
 
-function stringToTuple(T::DataType, s::String)
-    return Tuple(stringToVector(T,s))
+"""
+    get_stencils(parsed_toml, path...; centers)
+
+Same as [`read_stencils`](@ref)) but takes already parsed toml.
+"""
+get_stencils(parsed_toml, path...; centers) = get_stencils(parsed_toml[path[1]], path[2:end]...; centers=centers)
+function get_stencils(parsed_toml; centers)
+    @assert parsed_toml isa Vector{Vector{String}}
+    @assert length(centers) == length(parsed_toml)
+
+    stencils = ()
+    for i ∈ 1:length(parsed_toml)
+        stencil = get_stencil(parsed_toml[i], center = centers[i])
+        stencils = (stencils..., stencil)
+    end
+
+    return stencils
 end
 
-function stringToVector(T::DataType, s::String)
-    return T.(eval.(Meta.parse.(split(s))))
+"""
+    get_tuple(parsed_toml, path...)
+
+Same as [`read_tuple`](@ref)) but takes already parsed toml.
+"""
+get_tuple(parsed_toml, path...) = get_tuple(parsed_toml[path[1]], path[2:end]...)
+function get_tuple(parsed_toml)
+    @assert parsed_toml isa Vector{String}
+    t = Tuple(Float64.(parse_rational.(parsed_toml)))
+    return t
+end
+
+# TODO: Probably should be deleted once we have gotten rid of read_D2_operator()
+function toml_string_array_to_tuple(::Type{T}, arr::AbstractVector{String}) where T
+    return Tuple(T.(parse_rational.(arr)))
+end
+
+function parse_rational(str)
+    expr = Meta.parse(replace(str, "/"=>"//"))
+    return eval(:(Rational($expr)))
 end
 
 function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}

--- a/src/SbpOperators/stencil.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/src/SbpOperators/stencil.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -8,6 +8,29 @@
     end
 end
 
+"""
+    Stencil(weights::NTuple; center::Int)
+
+Create a stencil with the given weights with element `center` as the center of the stencil.
+"""
+function Stencil(weights::NTuple; center::Int)
+    N = length(weights)
+    range = (1, N) .- center
+
+    return Stencil(range, weights)
+end
+
+"""
+    scale(s::Stencil, a)
+
+Scale the weights of the stencil `s` with `a` and return a new stencil.
+"""
+function scale(s::Stencil, a)
+    return Stencil(s.range, a.*s.weights)
+end
+
+Base.eltype(::Stencil{T}) where T = T
+
 function flip(s::Stencil)
     range = (-s.range[2], -s.range[1])
     return Stencil(range, reverse(s.weights))
@@ -16,7 +39,7 @@
 # Provides index into the Stencil based on offset for the root element
 @inline function Base.getindex(s::Stencil, i::Int)
     @boundscheck if i < s.range[1] || s.range[2] < i
-        return eltype(s.weights)(0)
+        return zero(eltype(s))
     end
     return s.weights[1 + i - s.range[1]]
 end

--- a/test/Manifest.toml	Mon Nov 30 16:28:32 2020 +0100
+++ b/test/Manifest.toml	Fri Jan 01 16:34:55 2021 +0100
@@ -109,6 +109,12 @@
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "0.10.3"
 
+[[TOML]]
+deps = ["Dates"]
+git-tree-sha1 = "d0ac7eaad0fb9f6ba023a1d743edca974ae637c4"
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
+
 [[Test]]
 deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

--- a/test/Project.toml	Mon Nov 30 16:28:32 2020 +0100
+++ b/test/Project.toml	Fri Jan 01 16:34:55 2021 +0100
@@ -1,5 +1,6 @@
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

--- a/test/testDiffOps.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/test/testDiffOps.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -8,7 +8,7 @@
 @testset "DiffOps" begin
 #
 # @testset "BoundaryValue" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 #     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
 #
 #     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
@@ -69,7 +69,7 @@
 # end
 #
 # @testset "NormalDerivative" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 #     g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
 #
 #     d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
@@ -146,7 +146,7 @@
 # end
 #
 # @testset "BoundaryQuadrature" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 #     g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
 #
 #     H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())

--- a/test/testSbpOperators.jl	Mon Nov 30 16:28:32 2020 +0100
+++ b/test/testSbpOperators.jl	Fri Jan 01 16:34:55 2021 +0100
@@ -4,11 +4,112 @@
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 using LinearAlgebra
+using TOML
+
+import Sbplib.SbpOperators.Stencil
 
 @testset "SbpOperators" begin
 
+@testset "Stencil" begin
+    s = Stencil((-2,2), (1.,2.,2.,3.,4.))
+    @test s isa Stencil{Float64, 5}
+
+    @test eltype(s) == Float64
+    @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
+
+    @test Stencil((1,2,3,4), center=1) == Stencil((0, 3),(1,2,3,4))
+    @test Stencil((1,2,3,4), center=2) == Stencil((-1, 2),(1,2,3,4))
+    @test Stencil((1,2,3,4), center=4) == Stencil((-3, 0),(1,2,3,4))
+end
+
+@testset "parse_rational" begin
+    @test SbpOperators.parse_rational("1") isa Rational
+    @test SbpOperators.parse_rational("1") == 1//1
+    @test SbpOperators.parse_rational("1/2") isa Rational
+    @test SbpOperators.parse_rational("1/2") == 1//2
+    @test SbpOperators.parse_rational("37/13") isa Rational
+    @test SbpOperators.parse_rational("37/13") == 37//13
+end
+
+@testset "readoperator" begin
+    toml_str = """
+        [meta]
+        type = "equidistant"
+
+        [order2]
+        H.inner = ["1"]
+
+        D1.inner_stencil = ["-1/2", "0", "1/2"]
+        D1.closure_stencils = [
+            ["-1", "1"],
+        ]
+
+        d1.closure = ["-3/2", "2", "-1/2"]
+
+        [order4]
+        H.closure = ["17/48", "59/48", "43/48", "49/48"]
+
+        D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
+        D2.closure_stencils = [
+            [     "2",    "-5",      "4",       "-1",     "0",     "0"],
+            [     "1",    "-2",      "1",        "0",     "0",     "0"],
+            [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"],
+            [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"],
+        ]
+    """
+
+    parsed_toml = TOML.parse(toml_str)
+    @testset "get_stencil" begin
+        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil((-1/2, 0., 1/2), center=2)
+        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil((-1/2, 0., 1/2); center=1)
+        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil((-1/2, 0., 1/2); center=3)
+
+        @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil((1.,), center=1)
+
+        @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
+        @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
+    end
+
+    @testset "get_stencils" begin
+        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil((-1., 1.), center=1),)
+        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil((-1., 1.), center=2),)
+        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil((-1., 1.), center=2),)
+
+        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
+            Stencil((    2.,    -5.,      4.,     -1.,    0.,    0.), center=1),
+            Stencil((    1.,    -2.,      1.,      0.,    0.,    0.), center=1),
+            Stencil(( -4/43,  59/43, -110/43,   59/43, -4/43,    0.), center=1),
+            Stencil(( -1/49,     0.,   59/49, -118/49, 64/49, -4/49), center=1),
+        )
+
+        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
+            Stencil((    2.,    -5.,      4.,     -1.,    0.,    0.), center=4),
+            Stencil((    1.,    -2.,      1.,      0.,    0.,    0.), center=2),
+            Stencil(( -4/43,  59/43, -110/43,   59/43, -4/43,    0.), center=3),
+            Stencil(( -1/49,     0.,   59/49, -118/49, 64/49, -4/49), center=1),
+        )
+
+        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
+            Stencil((    2.,    -5.,      4.,     -1.,    0.,    0.), center=1),
+            Stencil((    1.,    -2.,      1.,      0.,    0.,    0.), center=2),
+            Stencil(( -4/43,  59/43, -110/43,   59/43, -4/43,    0.), center=3),
+            Stencil(( -1/49,     0.,   59/49, -118/49, 64/49, -4/49), center=4),
+        )
+
+        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
+        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4))
+        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2))
+    end
+
+    @testset "get_tuple" begin
+        @test get_tuple(parsed_toml, "order2", "d1", "closure") == (-3/2, 2, -1/2)
+
+        @test_throws AssertionError get_tuple(parsed_toml, "meta", "type")
+    end
+end
+
 # @testset "apply_quadrature" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 #     h = 0.5
 #
 #     @test apply_quadrature(op, h, 1.0, 10, 100) == h
@@ -29,7 +130,7 @@
 # end
 
 @testset "SecondDerivative" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     L = 3.5
     g = EquidistantGrid(101, 0.0, L)
     Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils)
@@ -69,7 +170,7 @@
 
 
 @testset "Laplace2D" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     Lx = 1.5
     Ly = 3.2
     g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
@@ -111,7 +212,7 @@
 end
 
 @testset "DiagonalQuadrature" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     Lx = π/2.
     Ly = Float64(π)
     g_1D = EquidistantGrid(77, 0.0, Lx)
@@ -169,14 +270,14 @@
         end
         # TODO: Bug in readOperator for 2nd order
         # # 2nd order
-        # op2 = readOperator(sbp_operators_path()*"d2_2nd.txt",sbp_operators_path()*"h_2nd.txt")
+        # op2 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
         # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure)
         # for i = 1:3
         #     @test integral(H2,v[i]) ≈ v[i+1] rtol = 1e-14
         # end
 
         # 4th order
-        op4 = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+        op4 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         H4 = diagonal_quadrature(g_1D,op4.quadratureClosure)
         for i = 1:4
             @test integral(H4,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
@@ -196,7 +297,7 @@
 end
 
 @testset "InverseDiagonalQuadrature" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     Lx = π/2.
     Ly = Float64(π)
     g_1D = EquidistantGrid(77, 0.0, Lx)
@@ -257,70 +358,147 @@
         @inferred Hi_xy'*v_2D
     end
 end
-#
-# @testset "BoundaryValue" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-#     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
-#
-#     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
-#     e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
-#     e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
-#     e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = zeros(Float64, 4, 5)
-#     v[:,5] = [1, 2, 3,4]
-#     v[:,4] = [1, 2, 3,4]
-#     v[:,3] = [4, 5, 6, 7]
-#     v[:,2] = [7, 8, 9, 10]
-#     v[:,1] = [10, 11, 12, 13]
-#
-#     @test e_w  isa TensorMapping{T,2,1} where T
-#     @test e_w' isa TensorMapping{T,1,2} where T
-#
-#     @test domain_size(e_w, (3,2)) == (2,)
-#     @test domain_size(e_e, (3,2)) == (2,)
-#     @test domain_size(e_s, (3,2)) == (3,)
-#     @test domain_size(e_n, (3,2)) == (3,)
-#
-#     @test size(e_w'*v) == (5,)
-#     @test size(e_e'*v) == (5,)
-#     @test size(e_s'*v) == (4,)
-#     @test size(e_n'*v) == (4,)
-#
-#     @test e_w'*v == [10,7,4,1.0,1]
-#     @test e_e'*v == [13,10,7,4,4.0]
-#     @test e_s'*v == [10,11,12,13.0]
-#     @test e_n'*v == [1,2,3,4.0]
-#
-#     g_x = [1,2,3,4.0]
-#     g_y = [5,4,3,2,1.0]
-#
-#     G_w = zeros(Float64, (4,5))
-#     G_w[1,:] = g_y
-#
-#     G_e = zeros(Float64, (4,5))
-#     G_e[4,:] = g_y
-#
-#     G_s = zeros(Float64, (4,5))
-#     G_s[:,1] = g_x
-#
-#     G_n = zeros(Float64, (4,5))
-#     G_n[:,5] = g_x
-#
-#     @test size(e_w*g_y) == (UnknownDim,5)
-#     @test size(e_e*g_y) == (UnknownDim,5)
-#     @test size(e_s*g_x) == (4,UnknownDim)
-#     @test size(e_n*g_x) == (4,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken e_w*g_y == G_w
-#     @test_broken e_e*g_y == G_e
-#     @test_broken e_s*g_x == G_s
-#     @test_broken e_n*g_x == G_n
-# end
+
+@testset "BoundaryRestrictrion" begin
+    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    g_1D = EquidistantGrid(11, 0.0, 1.0)
+    g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
+
+    @testset "Constructors" begin
+        @testset "1D" begin
+            e_l = BoundaryRestriction{Lower}(op.eClosure,size(g_1D)[1])
+            @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower())
+            @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
+            @test e_l isa TensorMapping{T,0,1} where T
+
+            e_r = BoundaryRestriction{Upper}(op.eClosure,size(g_1D)[1])
+            @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper())
+            @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
+            @test e_r isa TensorMapping{T,0,1} where T
+        end
+
+        @testset "2D" begin
+            e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            @test e_w isa InflatedTensorMapping
+            @test e_w isa TensorMapping{T,1,2} where T
+        end
+    end
+
+    e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
+    e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+
+    e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
+    e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
+    e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
+    e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
+
+    @testset "Sizes" begin
+        @testset "1D" begin
+            @test domain_size(e_l) == (11,)
+            @test domain_size(e_r) == (11,)
+
+            @test range_size(e_l) == ()
+            @test range_size(e_r) == ()
+        end
+
+        @testset "2D" begin
+            @test domain_size(e_w) == (11,15)
+            @test domain_size(e_e) == (11,15)
+            @test domain_size(e_s) == (11,15)
+            @test domain_size(e_n) == (11,15)
+
+            @test range_size(e_w) == (15,)
+            @test range_size(e_e) == (15,)
+            @test range_size(e_s) == (11,)
+            @test range_size(e_n) == (11,)
+        end
+    end
+
+
+    @testset "Application" begin
+        @testset "1D" begin
+            v = evalOn(g_1D,x->1+x^2)
+            u = fill(3.124)
+            @test (e_l*v)[] == v[1]
+            @test (e_r*v)[] == v[end]
+            @test (e_r*v)[1] == v[end]
+            @test e_l'*u == [u[]; zeros(10)]
+            @test e_r'*u == [zeros(10); u[]]
+        end
+
+        @testset "2D" begin
+            v = rand(11, 15)
+            u = fill(3.124)
+
+            @test e_w*v == v[1,:]
+            @test e_e*v == v[end,:]
+            @test e_s*v == v[:,1]
+            @test e_n*v == v[:,end]
+
+
+           g_x = rand(11)
+           g_y = rand(15)
+
+           G_w = zeros(Float64, (11,15))
+           G_w[1,:] = g_y
+
+           G_e = zeros(Float64, (11,15))
+           G_e[end,:] = g_y
+
+           G_s = zeros(Float64, (11,15))
+           G_s[:,1] = g_x
+
+           G_n = zeros(Float64, (11,15))
+           G_n[:,end] = g_x
+
+           @test e_w'*g_y == G_w
+           @test e_e'*g_y == G_e
+           @test e_s'*g_x == G_s
+           @test e_n'*g_x == G_n
+       end
+
+       @testset "Regions" begin
+           u = fill(3.124)
+           @test (e_l'*u)[Index(1,Lower)] == 3.124
+           @test (e_l'*u)[Index(2,Lower)] == 0
+           @test (e_l'*u)[Index(6,Interior)] == 0
+           @test (e_l'*u)[Index(10,Upper)] == 0
+           @test (e_l'*u)[Index(11,Upper)] == 0
+
+           @test (e_r'*u)[Index(1,Lower)] == 0
+           @test (e_r'*u)[Index(2,Lower)] == 0
+           @test (e_r'*u)[Index(6,Interior)] == 0
+           @test (e_r'*u)[Index(10,Upper)] == 0
+           @test (e_r'*u)[Index(11,Upper)] == 3.124
+       end
+    end
+
+    @testset "Inferred" begin
+        v = ones(Float64, 11)
+        u = fill(1.)
+
+        @inferred apply(e_l, v)
+        @inferred apply(e_r, v)
+
+        @inferred apply_transpose(e_l, u, 4)
+        @inferred apply_transpose(e_l, u, Index(1,Lower))
+        @inferred apply_transpose(e_l, u, Index(2,Lower))
+        @inferred apply_transpose(e_l, u, Index(6,Interior))
+        @inferred apply_transpose(e_l, u, Index(10,Upper))
+        @inferred apply_transpose(e_l, u, Index(11,Upper))
+
+        @inferred apply_transpose(e_r, u, 4)
+        @inferred apply_transpose(e_r, u, Index(1,Lower))
+        @inferred apply_transpose(e_r, u, Index(2,Lower))
+        @inferred apply_transpose(e_r, u, Index(6,Interior))
+        @inferred apply_transpose(e_r, u, Index(10,Upper))
+        @inferred apply_transpose(e_r, u, Index(11,Upper))
+    end
+
+end
 #
 # @testset "NormalDerivative" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 #     g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
 #
 #     d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
@@ -397,7 +575,7 @@
 # end
 #
 # @testset "BoundaryQuadrature" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 #     g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
 #
 #     H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())