--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Fri Jul 02 14:23:33 2021 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jan 19 14:44:24 2022 +0100
@@ -1,16 +1,16 @@
 """
     Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim}
-    Laplace(grid::AbstractGrid, fn; order)
+    Laplace(grid, filename; order)
 
 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
 `TensorMapping`. Additionally, `Laplace` stores the inner product and boundary
-operators relevant for constructing a SBP finite difference scheme as `TensorMapping`s.
+operators relevant for constructing a SBP finite difference scheme as a `TensorMapping`.
 
-Laplace(grid, fn; order) creates the Laplace operator defined on `grid`,
+`Laplace(grid, filename; order)` creates the Laplace operator defined on `grid`,
 where the operators are read from TOML. The differential operator is created
-using `laplace(grid::AbstractGrid,...)`.
+using `laplace(grid,...)`.
 
-Note that all properties of Laplace, excluding the Differential operator `D`, are
+Note that all properties of Laplace, excluding the differential operator `Laplace.D`, are
 abstract types. For performance reasons, they should therefore be
 accessed via the provided getter functions (e.g `inner_product(::Laplace)`).
 
@@ -21,32 +21,34 @@
     H_inv::TensorMapping # Inverse inner product operator
     e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators.
     d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators
-    H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product?
+    H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators
 end
 export Laplace
 
-function Laplace(grid::AbstractGrid, fn; order)
-    # TODO: Removed once we can construct the volume and
-    # boundary operators by op(grid, fn; order,...).
+function Laplace(grid, filename; order)
+    
     # Read stencils
-    op = read_D2_operator(fn; order)
-    D_inner_stecil = op.innerStencil
-    D_closure_stencils = op.closureStencils
-    H_closure_stencils = op.quadratureClosure
-    e_closure_stencil = op.eClosure
-    d_closure_stencil = op.dClosure
+    stencil_set = read_stencil_set(filename; order)
+    # TODO: Removed once we can construct the volume and
+    # boundary operators by op(grid, read_stencil_set(fn; order,...)).
+    D_inner_stecil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+    D_closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+    H_inner_stencils = parse_scalar(stencil_set["H"]["inner"])
+    H_closure_stencils = parse_tuple(stencil_set["H"]["closure"])
+    e_closure_stencil = parse_stencil(stencil_set["e"]["closure"])
+    d_closure_stencil = parse_stencil(stencil_set["d1"]["closure"])
 
     # Volume operators
     Δ =  laplace(grid, D_inner_stecil, D_closure_stencils)
-    H =  inner_product(grid, H_closure_stencils)
-    H⁻¹ =  inverse_inner_product(grid, H_closure_stencils)
+    H =  inner_product(grid, H_inner_stencils, H_closure_stencils)
+    H⁻¹ = inverse_inner_product(grid, H_inner_stencils, H_closure_stencils)
 
     # Boundary operator - id pairs
     ids = boundary_identifiers(grid)
     n_ids = length(ids)
-    e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid,e_closure_stencil,ids[i]),n_ids)
-    d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids)
-    Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids)
+    e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid, e_closure_stencil, ids[i]), n_ids)
+    d_pairs = ntuple(i -> ids[i] => normal_derivative(grid, d_closure_stencil, ids[i]), n_ids)
+    Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid, ids[i]), H_inner_stencils, H_closure_stencils), n_ids)
 
     return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs))
 end
@@ -60,7 +62,7 @@
 
 
 """
-    inner_product(L::Lapalace)
+    inner_product(L::Laplace)
 
 Returns the inner product operator associated with `L`
 
@@ -70,7 +72,7 @@
 
 
 """
-    inverse_inner_product(L::Lapalace)
+    inverse_inner_product(L::Laplace)
 
 Returns the inverse of the inner product operator associated with `L`
 
@@ -80,65 +82,64 @@
 
 
 """
-    boundary_restriction(L::Lapalace,id::BoundaryIdentifier)
-    boundary_restriction(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
-    boundary_restriction(L::Lapalace,ids...)
+    boundary_restriction(L::Laplace, id::BoundaryIdentifier)
+    boundary_restriction(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
+    boundary_restriction(L::Laplace, ids...)
 
 Returns boundary restriction operator(s) associated with `L` for the boundary(s)
 identified by id(s).
 
 """
-boundary_restriction(L::Laplace,id::BoundaryIdentifier) = L.e[id]
-boundary_restriction(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N)
-boundary_restriction(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N)
+boundary_restriction(L::Laplace, id::BoundaryIdentifier) = L.e[id]
+boundary_restriction(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N)
+boundary_restriction(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N)
 export boundary_restriction
 
 
 """
-    normal_derivative(L::Lapalace,id::BoundaryIdentifier)
-    normal_derivative(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
-    normal_derivative(L::Lapalace,ids...)
+    normal_derivative(L::Laplace, id::BoundaryIdentifier)
+    normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
+    normal_derivative(L::Laplace, ids...)
 
 Returns normal derivative operator(s) associated with `L` for the boundary(s)
 identified by id(s).
 
 """
-normal_derivative(L::Laplace,id::BoundaryIdentifier) = L.d[id]
-normal_derivative(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
-normal_derivative(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
+normal_derivative(L::Laplace, id::BoundaryIdentifier) = L.d[id]
+normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
+normal_derivative(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
 export normal_derivative
 
 
-# TODO: boundary_inner_product?
 """
-    boundary_quadrature(L::Lapalace,id::BoundaryIdentifier)
-    boundary_quadrature(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
-    boundary_quadrature(L::Lapalace,ids...)
+    boundary_quadrature(L::Laplace, id::BoundaryIdentifier)
+    boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
+    boundary_quadrature(L::Laplace, ids...)
 
 Returns boundary quadrature operator(s) associated with `L` for the boundary(s)
 identified by id(s).
 
 """
-boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id]
-boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
-boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
+boundary_quadrature(L::Laplace, id::BoundaryIdentifier) = L.H_boundary[id]
+boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
+boundary_quadrature(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
 export boundary_quadrature
 
 
 """
-    laplace(grid, inner_stencil, closure_stencils)
+    laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
 
 Creates the Laplace operator operator `Δ` as a `TensorMapping`
 
-`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,N on the N-dimensional
-`grid`, using the stencil `inner_stencil` in the interior and a set of stencils
-`closure_stencils` for the points in the closure regions.
+`Δ` 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.
 
 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.
 """
-function laplace(grid::AbstractGrid, inner_stencil, closure_stencils)
+function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
     Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
     for d = 2:dimension(grid)
         Δ += second_derivative(grid, inner_stencil, closure_stencils, d)