changeset 1291:356ec6a72974 refactor/grids

Implement changes in SbpOperators
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Mar 2023 09:48:00 +0100
parents 31d0b7638304
children 6753b210d0ab
files src/SbpOperators/volumeops/constant_interior_scaling_operator.jl src/SbpOperators/volumeops/derivatives/dissipation.jl src/SbpOperators/volumeops/derivatives/first_derivative.jl src/SbpOperators/volumeops/derivatives/second_derivative.jl src/SbpOperators/volumeops/inner_products/inner_product.jl src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl src/SbpOperators/volumeops/laplace/laplace.jl src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl src/SbpOperators/volumeops/volume_operator.jl test/Grids/grid_test.jl test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl test/SbpOperators/volumeops/derivatives/dissipation_test.jl test/SbpOperators/volumeops/derivatives/first_derivative_test.jl test/SbpOperators/volumeops/derivatives/second_derivative_test.jl test/SbpOperators/volumeops/inner_products/inner_product_test.jl test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl test/SbpOperators/volumeops/laplace/laplace_test.jl test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl test/SbpOperators/volumeops/volume_operator_test.jl
diffstat 19 files changed, 220 insertions(+), 362 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Tue Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -18,7 +18,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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -10,30 +10,31 @@
 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(g::TensorGrid, p, direction)
+    op = undivided_skewed04(g.grids[direction], p)
+    return LazyTensors.inflate(op, 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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -1,48 +1,36 @@
 """
-    first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+    first_derivative(g::TensorGrid, stencil_set, direction)
 
 Creates the first-derivative operator `D1` as a `LazyTensor`
 
-`D1` approximates the first-derivative d/dξ on `grid` along the coordinate dimension specified by
+`D1` approximates the first-derivative d/dξ on `g` 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.
 
-On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the inflation of
-a `VolumeOperator`.
-
 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)
 
-Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`.
+Creates a `first_derivative` operator on a `EquidistantGrid` given a `StencilSet`.
 """
-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, closure_stencils)
 
-Creates a `first_derivative` operator on a 1D `grid` given a `stencil_set`.
+Creates a `first_derivative` operator on a `EquidistantGrid` given an `inner_stencil` and a`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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -1,48 +1,36 @@
 """
-    second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+    second_derivative(g::EquidistantGrid, inner_stencil, closure_stencils, direction)
 
 Creates the second-derivative operator `D2` as a `LazyTensor`
 
-`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by
+`D2` approximates the second-derivative d²/dξ² on `g` 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.
 
-On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the inflation of
-a `VolumeOperator`.
-
 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, stencil_set)
 
-Creates a `second_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`.
+Creates a `second_derivative` operator on a 1D `g` given a `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, inner_stencil, closure_stencils)
 
-Creates a `second_derivative` operator on a 1D `grid` given a `stencil_set`.
+Creates a `second_derivative` operator on a 1D `g` 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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -15,32 +15,17 @@
 
 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(tg::TensorGrid, stencil_set::StencilSet)
+    return mapreduce(g->inner_product(g,stencil_set), ⊗, tg.grids)
 end
 
-function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
-    h = spacing(grid)[1]
+function inner_product(g::EquidistantGrid, stencil_set::StencilSet)
+    interior_weight = parse_scalar(stencil_set["H"]["inner"])
+    closure_weights = parse_tuple(stencil_set["H"]["closure"])
 
-    H = SbpOperators.ConstantInteriorScalingOperator(grid, h*interior_weight, h.*closure_weights)
-    return H
+    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) = IdentityTensor{component_type(g)}()
 
-Creates a `inner_product` operator on `grid` given a `stencil_set`.
-"""
-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
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -12,31 +12,16 @@
 
 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(tg::TensorGrid, stencil_set::StencilSet)
+    return mapreduce(g->inverse_inner_product(g,stencil_set), ⊗, tg.grids)
 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⁻¹
+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"])
+
+    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)
-
-Creates a `inverse_inner_product` operator on `grid` given a `stencil_set`.
-"""
-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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -3,7 +3,7 @@
 
 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
 `LazyTensor`. Additionally `Laplace` stores the `StencilSet`
-used to construct the `LazyTensor `.
+used to construct the `LazyTensor`.
 """
 struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim}
     D::TM       # Difference operator
@@ -17,11 +17,9 @@
 
 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 +30,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 the given grid
 
-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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -1,22 +1,3 @@
-"""
-    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}
 
@@ -37,7 +18,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 Mar 07 09:21:27 2023 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -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/grid_test.jl	Tue Mar 07 09:21:27 2023 +0100
+++ b/test/Grids/grid_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -52,6 +52,8 @@
     # Multi-argument functions
     f(x,y) = sin(x)*cos(y)
     @test eval_on(g, f) == map(x̄->f(x̄...), g)
+
+    #TODO: inference test!
 end
 
 @testset "_ncomponents" begin
--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Tue Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -5,7 +5,7 @@
 import Sbplib.SbpOperators: ConstantInteriorScalingOperator
 using Sbplib.Grids
 
-@test_skip @testset "ConstantInteriorScalingOperator" begin
+@testset "ConstantInteriorScalingOperator" begin
     @test ConstantInteriorScalingOperator(1, (2,3), 10) isa ConstantInteriorScalingOperator{Int,2}
     @test ConstantInteriorScalingOperator(1., (2.,3.), 10) isa ConstantInteriorScalingOperator{Float64,2}
 
@@ -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 Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -26,8 +26,8 @@
     x^k/factorial(k)
 end
 
-@test_skip @testset "undivided_skewed04" begin
-    g = EquidistantGrid(20, 0., 11.)
+@testset "undivided_skewed04" begin
+    g = equidistant_grid(20, 0., 11.)
     D,Dᵀ = undivided_skewed04(g, 1)
 
     @test D isa LazyTensor{Float64,1,1}
@@ -35,18 +35,20 @@
 
      @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
         end
+
+        # TODO: Add 2D tests
     end
 
     @testset "transpose equality" begin
@@ -67,7 +69,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)
 
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -21,76 +21,74 @@
 end
 
 @testset "first_derivative" begin
-    @test_skip @testset "Constructors" begin
+    @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
 
-    @test_skip @testset "Accuracy conditions" begin
+    @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
 
-    @test_skip @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 "Accuracy on function" begin
+        @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 Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -8,31 +8,25 @@
 
 # TODO: Refactor these test to look more like the tests in first_derivative_test.jl.
 
-@test_skip @testset "SecondDerivative" begin
+@testset "SecondDerivative" begin
     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"])
     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 Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -6,51 +6,43 @@
 
 import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
-@test_skip @testset "Diagonal-stencil inner_product" begin
+@testset "Diagonal-stencil inner_product" begin
     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 Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -6,32 +6,25 @@
 
 import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
-@test_skip @testset "Diagonal-stencil inverse_inner_product" begin
+@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 Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -4,40 +4,40 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
-@test_skip @testset "Laplace" begin
+@testset "Laplace" begin
     # 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.))
+    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 "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.
@@ -66,27 +66,25 @@
     end
 end
 
-@test_skip @testset "laplace" begin
+@testset "laplace" begin
     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.))
+    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 Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -6,40 +6,9 @@
 
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.StencilOperatorDistinctClosures
-import Sbplib.SbpOperators.stencil_operator_distinct_closures
 
-@test_skip @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
-
-@test_skip @testset "StencilOperatorDistinctClosures" begin
-    g = EquidistantGrid(11, 0., 1.)
+@testset "StencilOperatorDistinctClosures" begin
+    g = equidistant_grid(11, 0., 1.)
 
     lower_closure = (
         Stencil(-1,1, center=1),
--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Tue Mar 07 09:21:27 2023 +0100
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Tue Mar 07 09:48:00 2023 +0100
@@ -11,10 +11,10 @@
 import Sbplib.SbpOperators.even
 
 
-@test_skip @testset "VolumeOperator" begin
+@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)