Mercurial > repos > public > sbplib_julia

changeset 702:3cd582257072 feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 15 Feb 2021 11:30:34 +0100
parents 54ce3f9771e5 (current diff) 38f9894279cd (diff)
children 988e9cfcd58d
files Notes.md src/SbpOperators/volumeops/laplace/laplace.jl src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl src/SbpOperators/volumeops/quadratures/quadrature.jl test/testSbpOperators.jl
diffstat 15 files changed, 353 insertions(+), 314 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Sat Feb 06 15:26:14 2021 +0100
+++ b/Notes.md	Mon Feb 15 11:30:34 2021 +0100
@@ -299,3 +299,16 @@
 component(gf,:,2) # Andra kolumnen av en matris
 @ourview gf[:,:][2]
 ```
+
+## Grids embedded in higher dimensions
+
+For grids generated by asking for boundary grids for a regular grid, it would
+make sense if these grids knew they were embedded in a higher dimension. They
+would return coordinates in the full room. This would make sense when
+drawing points for example, or when evaluating functions on the boundary.
+
+Implementation of this is an issue that requires some thought. Adding an extra
+"Embedded" type for each grid would make it easy to understand each type but
+contribute to "type bloat". On the other hand adapting existing types to
+handle embeddedness would complicate the now very simple grid types. Are there
+other ways of doing the implentation?
--- a/src/Grids/EquidistantGrid.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/Grids/EquidistantGrid.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -1,11 +1,19 @@
 """
-    EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}
+    EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T})
+	EquidistantGrid{T}()
+
+`EquidistantGrid` is a grid with equidistant grid spacing per coordinat direction.
 
-EquidistantGrid is a grid with equidistant grid spacing per coordinat direction.
-The domain is defined through the two points P1 = x̄₁, P2 = x̄₂ by the exterior
-product of the vectors obtained by projecting (x̄₂-x̄₁) onto the coordinate
-directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2) the domain is defined
-as (-1,1)x(0,2). The side lengths of the grid are not allowed to be negative
+`EquidistantGrid(size, limit_lower, limit_upper)` construct the grid with the
+domain defined by the two points P1, and P2 given by `limit_lower` and
+`limit_upper`. The length of the domain sides are given by the components of
+(P2-P1). E.g for a 2D grid with P1=(-1,0) and P2=(1,2) the domain is defined
+as (-1,1)x(0,2). The side lengths of the grid are not allowed to be negative.
+The number of equidistantly spaced points in each coordinate direction are given
+by `size`.
+
+`EquidistantGrid{T}()` constructs a 0-dimensional grid.
+
 """
 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
     size::NTuple{Dim, Int}
@@ -22,6 +30,9 @@
         end
         return new{Dim,T}(size, limit_lower, limit_upper)
     end
+
+	# Specialized constructor for 0-dimensional grid
+	EquidistantGrid{T}() where T = new{0,T}((),(),())
 end
 export EquidistantGrid
 
@@ -35,9 +46,9 @@
 	return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
 end
 
-function Base.eachindex(grid::EquidistantGrid)
-    CartesianIndices(grid.size)
-end
+Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
+
+Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
 
 Base.size(g::EquidistantGrid) = g.size
 
@@ -104,3 +115,23 @@
 """
 boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),dimension(g)))...)...,)
 export boundary_identifiers
+
+
+"""
+    boundary_grid(grid::EquidistantGrid,id::CartesianBoundary)
+	boundary_grid(::EquidistantGrid{1},::CartesianBoundary{1})
+
+Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
+orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
+grid is a zero-dimensional grid.
+"""
+function boundary_grid(grid::EquidistantGrid,id::CartesianBoundary)
+	dims = collect(1:dimension(grid))
+	orth_dims = dims[dims .!= dim(id)]
+	if orth_dims == dims
+		throw(DomainError("boundary identifier not matching grid"))
+	end
+    return restrict(grid,orth_dims)
+end
+export boundary_grid
+boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
--- a/src/SbpOperators/SbpOperators.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -10,8 +10,8 @@
 include("volumeops/volume_operator.jl")
 include("volumeops/derivatives/secondderivative.jl")
 include("volumeops/laplace/laplace.jl")
-include("volumeops/quadratures/quadrature.jl")
-include("volumeops/quadratures/inverse_quadrature.jl")
+include("volumeops/inner_products/inner_product.jl")
+include("volumeops/inner_products/inverse_inner_product.jl")
 include("boundaryops/boundary_operator.jl")
 include("boundaryops/boundary_restriction.jl")
 include("boundaryops/normal_derivative.jl")
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -1,6 +1,6 @@
 """
-    BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
 
 Creates the boundary restriction operator `e` as a `TensorMapping`
 
@@ -9,7 +9,7 @@
 On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
 a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
 """
-BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
-BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = BoundaryRestriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
+boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
 
-export BoundaryRestriction
+export boundary_restriction
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -1,6 +1,6 @@
 """
-    NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
 
 Creates the normal derivative boundary operator `d` as a `TensorMapping`
 
@@ -9,10 +9,10 @@
 On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
 a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
 """
-function NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+function normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
 end
-NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = NormalDerivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
-export NormalDerivative
+normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+export normal_derivative
--- a/src/SbpOperators/readoperator.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/readoperator.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -26,13 +26,13 @@
         closureStencils = (closureStencils..., get_stencil(operators, "D2", "closure_stencils", i; center=i))
     end
     # TODO: Get rid of the padding here. Any padding should be handled by the consturctor accepting the stencils.
-    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)
+    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)
 
     q_tuple = pad_tuple(toml_string_array_to_tuple(Float64, H["closure"]), boundarySize)
     quadratureClosure = Vector{typeof(innerStencil)}()
     for i ∈ 1:boundarySize
-        quadratureClosure = (quadratureClosure..., Stencil((q_tuple[i],), center=1))
+        quadratureClosure = (quadratureClosure..., Stencil(q_tuple[i], center=1))
     end
 
     d2 = SbpOperators.D2(
@@ -99,7 +99,7 @@
         center = div(width,2)+1
     end
 
-    return Stencil(Tuple(stencil_weights), center=center)
+    return Stencil(stencil_weights..., center=center)
 end
 
 """
--- a/src/SbpOperators/stencil.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/stencil.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -1,3 +1,5 @@
+export CenteredStencil
+
 struct Stencil{T<:Real,N}
     range::Tuple{Int,Int}
     weights::NTuple{N,T}
@@ -13,13 +15,24 @@
 
 Create a stencil with the given weights with element `center` as the center of the stencil.
 """
-function Stencil(weights::NTuple; center::Int)
+function Stencil(weights::Vararg{Number}; center::Int)
     N = length(weights)
     range = (1, N) .- center
 
     return Stencil(range, weights)
 end
 
+function CenteredStencil(weights::Vararg)
+    if iseven(length(weights))
+        throw(ArgumentError("a centered stencil must have an odd number of weights."))
+    end
+
+    r = length(weights) ÷ 2
+
+    return Stencil((-r, r), weights)
+end
+
+
 """
     scale(s::Stencil, a)
 
--- a/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -1,6 +1,6 @@
 """
-    SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
-    SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
+    second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
+    second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
 
 Creates the second-derivative operator `D2` as a `TensorMapping`
 
@@ -12,9 +12,9 @@
 one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
 Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
 """
-function SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
+function second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
 end
-SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = SecondDerivative(grid,inner_stencil,closure_stencils,1)
-export SecondDerivative
+second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
+export second_derivative
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -0,0 +1,29 @@
+"""
+    inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
+
+Creates the discrete inner product operator `H` as a `TensorMapping` on an equidistant
+grid, defined as `(u,v)  = u'Hv` for grid functions `u,v`.
+
+`inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)` creates
+`H` on `grid` the using a set of stencils `closure_stencils` for the points in
+the closure regions and the stencil and `inner_stencil` in the interior. If
+`inner_stencil` is omitted a central interior stencil with weight 1 is used.
+
+On a 1-dimensional `grid`, `H` is a `VolumeOperator`. On a N-dimensional
+`grid`, `H` is the outer product of the 1-dimensional inner product operators in
+each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
+`H` is a 0-dimensional `IdentityMapping`.
+"""
+function inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil = CenteredStencil(one(eltype(grid))))
+    h = spacing(grid)
+    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
+    for i ∈ 2:dimension(grid)
+        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
+        H = H∘Hᵢ
+    end
+    return H
+end
+export inner_product
+
+inner_product(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -0,0 +1,43 @@
+"""
+    inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
+    inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})
+
+Creates the inverse inner product operator `H⁻¹` as a `TensorMapping` on an
+equidistant grid. `H⁻¹` is defined implicitly by `H⁻¹∘H = I`, where
+`H` is the corresponding inner product operator and `I` is the `IdentityMapping`.
+
+`inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)`
+constructs `H⁻¹` using a set of stencils `inv_closure_stencils` for the points
+in the closure regions and the stencil `inv_inner_stencil` in the interior. If
+`inv_closure_stencils` is omitted, a central interior stencil with weight 1 is used.
+
+`inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})`
+constructs a diagonal inverse inner product operator where `closure_stencils` are the
+closure stencils of `H` (not `H⁻¹`!).
+
+On a 1-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a N-dimensional
+`grid`, `H⁻¹` is the outer product of the 1-dimensional inverse inner product
+operators in each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
+`H⁻¹` is a 0-dimensional `IdentityMapping`.
+"""
+function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid))))
+    h⁻¹ = inverse_spacing(grid)
+    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
+    for i ∈ 2:dimension(grid)
+        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
+        H⁻¹ = H⁻¹∘Hᵢ⁻¹
+    end
+    return H⁻¹
+end
+export inverse_inner_product
+
+inverse_inner_product(grid::EquidistantGrid{0}, inv_closure_stencils, inv_inner_stencil) = IdentityMapping{eltype(grid)}()
+
+function inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T}
+     inv_closure_stencils = reciprocal_stencil.(closure_stencils)
+     inv_inner_stencil = CenteredStencil(one(T))
+     return inverse_inner_product(grid, inv_closure_stencils, inv_inner_stencil)
+end
+
+reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -37,14 +37,14 @@
 
     # Volume operators
     Δ =  laplace(grid, D_inner_stecil, D_closure_stencils)
-    H =  quadrature(grid, H_closure_stencils)
-    H⁻¹ =  InverseDiagonalQuadrature(grid, H_closure_stencils)
+    H =  inner_product(grid, H_closure_stencils)
+    H⁻¹ =  inverse_inner_product(grid, H_closure_stencils)
 
     # Boundary operator - id pairs
     bids = boundary_identifiers(grid)
-    e_pairs = ntuple(i -> Pair(bids[i],BoundaryRestriction(grid,e_closure_stencil,bids[i])),length(bids))
-    d_pairs = ntuple(i -> Pair(bids[i],NormalDerivative(grid,d_closure_stencil,bids[i])),length(bids))
-    Hᵧ_pairs = ntuple(i -> Pair(bids[i],boundary_quadrature(grid,H_closure_stencils,bids[i])),length(bids))
+    e_pairs = ntuple(i -> Pair(bids[i],boundary_restriction(grid,e_closure_stencil,bids[i])),length(bids))
+    d_pairs = ntuple(i -> Pair(bids[i],normal_derivative(grid,d_closure_stencil,bids[i])),length(bids))
+    Hᵧ_pairs = ntuple(i -> Pair(bids[i],inner_product(boundary_grid(grid,bids[i]),H_closure_stencils)),length(bids))
 
     return Laplace(Δ, H, H⁻¹, Dict(e_pairs), Dict(d_pairs), Dict(Hᵧ_pairs))
 end
@@ -73,13 +73,14 @@
 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 a `SecondDerivative`. On a multi-dimensional `grid`, `Δ` is the sum of
-multi-dimensional `SecondDerivative`s where the sum is carried out lazily.
+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::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
-    Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1)
+    Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
     for d = 2:Dim
-        Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d)
+        Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
     end
     return Δ
 end
--- a/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl	Sat Feb 06 15:26:14 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-
-"""
-    InverseQuadrature(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
-
-Creates the inverse `H⁻¹` of the quadrature operator as a `TensorMapping`
-
-The inverse quadrature approximates the integral operator on the grid using
-`inv_inner_stencil` in the interior and a set of stencils `inv_closure_stencils`
-for the points in the closure regions.
-
-On a one-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a multi-dimensional
-`grid`, `H` is the outer product of the 1-dimensional inverse quadrature operators in
-each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details.
-"""
-function InverseQuadrature(grid::EquidistantGrid{Dim}, inv_inner_stencil, inv_closure_stencils) where Dim
-    h⁻¹ = inverse_spacing(grid)
-    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
-    for i ∈ 2:Dim
-        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
-        H⁻¹ = H⁻¹∘Hᵢ⁻¹
-    end
-    return H⁻¹
-end
-export InverseQuadrature
-
-"""
-    InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
-
-Creates the inverse of the diagonal quadrature operator defined by the inner stencil
-1/h and a set of 1-element closure stencils in `closure_stencils`. Note that
-the closure stencils are those of the quadrature operator (and not the inverse).
-"""
-function InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {T,M}
-    inv_inner_stencil = Stencil(Tuple{T}(1),center=1)
-    inv_closure_stencils = reciprocal_stencil.(closure_stencils)
-    return InverseQuadrature(grid, inv_inner_stencil, inv_closure_stencils)
-end
-export InverseDiagonalQuadrature
-
-reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
--- a/src/SbpOperators/volumeops/quadratures/quadrature.jl	Sat Feb 06 15:26:14 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-"""
-    quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils)
-    quadrature(grid::EquidistantGrid, closure_stencils)
-
-Creates the quadrature operator `H` as a `TensorMapping`
-
-`H` approximiates the integral operator on `grid` the using the stencil
-`inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions. If `inner_stencil` is omitted a central
-interior stencil with weight 1 is used.
-
-On a one-dimensional `grid`, `H` is a `VolumeOperator`. On a multi-dimensional
-`grid`, `H` is the outer product of the 1-dimensional quadrature operators in
-each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details.
-"""
-function quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils) where Dim
-    h = spacing(grid)
-    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
-    for i ∈ 2:dimension(grid)
-        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
-        H = H∘Hᵢ
-    end
-    return H
-end
-export quadrature
-
-function quadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T}}) where {M,T}
-    inner_stencil = Stencil(Tuple{T}(1),center=1)
-    return quadrature(grid, inner_stencil, closure_stencils)
-end
-
-"""
-    boundary_quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils, id::CartesianBoundary)
-    boundary_quadrature(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, id)
-    boundary_quadrature(grid::EquidistantGrid, closure_stencils, id)
-
-Creates the lower-dimensional quadrature operator associated with the boundary
-of `grid` specified by `id`. The quadrature operator is defined on the grid
-spanned by the dimensions orthogonal to the boundary coordinate direction.
-If the dimension of `grid` is 1, then the boundary quadrature is the 0-dimensional
-`IdentityMapping`. If `inner_stencil` is omitted a central interior stencil with
-weight 1 is used.
-"""
-function boundary_quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils, id::CartesianBoundary)
-    return quadrature(orthogonal_grid(grid,dim(id)),inner_stencil,closure_stencils)
-end
-export boundary_quadrature
-
-function boundary_quadrature(grid::EquidistantGrid{1}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, id::CartesianBoundary{1}) where {M,T}
-    return IdentityMapping{T}()
-end
-
-function boundary_quadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T}}, id::CartesianBoundary) where {M,T}
-    inner_stencil = Stencil(Tuple{T}(1),center=1)
-    return boundary_quadrature(grid,inner_stencil,closure_stencils,id)
-end
-
-"""
-    orthogonal_grid(grid,dim)
-
-Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
-orthogonal to `dim`.
-"""
-function orthogonal_grid(grid,dim)
-    dims = collect(1:dimension(grid))
-    orth_dims = dims[dims .!= dim]
-    return restrict(grid,orth_dims)
-end
--- a/test/testGrids.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/test/testGrids.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -13,9 +13,12 @@
     @test_throws DomainError EquidistantGrid(1,1.0,-1.0)
     @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,))
 
-    # size
-    @test size(EquidistantGrid(4,0.0,1.0)) == (4,)
-    @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+    @testset "Base" begin
+        @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64
+        @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int
+        @test size(EquidistantGrid(4,0.0,1.0)) == (4,)
+        @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+    end
 
     # dimension
     @test dimension(EquidistantGrid(4,0.0,1.0)) == 1
@@ -63,6 +66,39 @@
         @test boundary_identifiers(g) == bids
         @inferred boundary_identifiers(g)
     end
+
+    @testset "boundary_grid" begin
+            @testset "1D" begin
+                g = EquidistantGrid(5,0.0,2.0)
+                (id_l, id_r) = boundary_identifiers(g)
+                @test boundary_grid(g,id_l) == EquidistantGrid{Float64}()
+                @test boundary_grid(g,id_r) == EquidistantGrid{Float64}()
+                @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}())
+                @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}())
+            end
+            @testset "2D" begin
+                g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
+                (id_w, id_e, id_s, id_n) = boundary_identifiers(g)
+                @test boundary_grid(g,id_w) == restrict(g,2)
+                @test boundary_grid(g,id_e) == restrict(g,2)
+                @test boundary_grid(g,id_s) == restrict(g,1)
+                @test boundary_grid(g,id_n) == restrict(g,1)
+                @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
+            end
+            @testset "3D" begin
+                g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
+                (id_w, id_e,
+                 id_s, id_n,
+                 id_t, id_b) = boundary_identifiers(g)
+                @test boundary_grid(g,id_w) == restrict(g,[2,3])
+                @test boundary_grid(g,id_e) == restrict(g,[2,3])
+                @test boundary_grid(g,id_s) == restrict(g,[1,3])
+                @test boundary_grid(g,id_n) == restrict(g,[1,3])
+                @test boundary_grid(g,id_t) == restrict(g,[1,2])
+                @test boundary_grid(g,id_b) == restrict(g,[1,2])
+                @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
+            end
+    end
 end
 
 end
--- a/test/testSbpOperators.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/test/testSbpOperators.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -24,9 +24,12 @@
     @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))
+    @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))
+
+    @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
+    @test_throws ArgumentError CenteredStencil(1,2,3,4)
 end
 
 @testset "parse_rational" begin
@@ -67,40 +70,40 @@
 
     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", "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 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, "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),
+            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),
+            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),
+            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))
@@ -116,8 +119,8 @@
 end
 
 @testset "VolumeOperator" begin
-    inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2)
-    closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2))
+    inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
+    closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2))
     g_1D = EquidistantGrid(11,0.,1.)
     g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
     g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
@@ -238,13 +241,13 @@
 
     @testset "Constructors" begin
         @testset "1D" begin
-            Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
-            @test Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1)
+            Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
+            @test Dₓₓ == second_derivative(g_1D,op.innerStencil,op.closureStencils,1)
             @test Dₓₓ isa VolumeOperator
         end
         @testset "2D" begin
-            Dₓₓ = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,1)
-            D2 = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+            Dₓₓ = second_derivative(g_2D,op.innerStencil,op.closureStencils,1)
+            D2 = second_derivative(g_1D,op.innerStencil,op.closureStencils)
             I = IdentityMapping{Float64}(size(g_2D)[2])
             @test Dₓₓ == D2⊗I
             @test Dₓₓ isa TensorMapping{T,2,2} where T
@@ -269,7 +272,7 @@
             # implies that L*v should be exact for monomials up to order 2.
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+                Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
                 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
                 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
                 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
@@ -280,7 +283,7 @@
             # implies that L*v should be exact for monomials up to order 3.
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+                Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
                 # NOTE: high tolerances for checking the "exact" differentiation
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
@@ -306,7 +309,7 @@
             # implies that L*v should be exact for binomials up to order 2.
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
+                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,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
@@ -317,7 +320,7 @@
             # implies that L*v should be exact for binomials up to order 3.
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
+                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
                 # NOTE: high tolerances for checking the "exact" differentiation
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
@@ -338,21 +341,21 @@
         @testset "1D" begin
             # Create all tensor mappings included in Laplace
             Δ = laplace(g_1D, op.innerStencil, op.closureStencils)
-            H = quadrature(g_1D, op.quadratureClosure)
-            Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure)
+            H = inner_product(g_1D, op.quadratureClosure)
+            Hi = inverse_inner_product(g_1D, op.quadratureClosure)
 
             (id_l, id_r) = boundary_identifiers(g_1D)
 
-            e_l = BoundaryRestriction(g_1D,op.eClosure,id_l)
-            e_r = BoundaryRestriction(g_1D,op.eClosure,id_r)
+            e_l = boundary_restriction(g_1D,op.eClosure,id_l)
+            e_r = boundary_restriction(g_1D,op.eClosure,id_r)
             e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r))
 
-            d_l = NormalDerivative(g_1D,op.dClosure,id_l)
-            d_r = NormalDerivative(g_1D,op.dClosure,id_r)
+            d_l = normal_derivative(g_1D,op.dClosure,id_l)
+            d_r = normal_derivative(g_1D,op.dClosure,id_r)
             d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r))
 
-            H_l = boundary_quadrature(g_1D,op.quadratureClosure,id_r)
-            H_r = boundary_quadrature(g_1D,op.quadratureClosure,id_r)
+            H_l = inner_product(boundary_grid(g_1D,id_l),op.quadratureClosure)
+            H_r = inner_product(boundary_grid(g_1D,id_r),op.quadratureClosure)
             Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r))
 
             # TODO: Not sure why this doesnt work? Comparing the fields of
@@ -372,37 +375,37 @@
         @testset "3D" begin
             # Create all tensor mappings included in Laplace
             Δ = laplace(g_3D, op.innerStencil, op.closureStencils)
-            H = quadrature(g_3D, op.quadratureClosure)
-            Hi = InverseDiagonalQuadrature(g_3D, op.quadratureClosure)
+            H = inner_product(g_3D, op.quadratureClosure)
+            Hi = inverse_inner_product(g_3D, op.quadratureClosure)
 
             (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
 
-            e_l = BoundaryRestriction(g_3D,op.eClosure,id_l)
-            e_r = BoundaryRestriction(g_3D,op.eClosure,id_r)
-            e_s = BoundaryRestriction(g_3D,op.eClosure,id_s)
-            e_n = BoundaryRestriction(g_3D,op.eClosure,id_n)
-            e_b = BoundaryRestriction(g_3D,op.eClosure,id_b)
-            e_t = BoundaryRestriction(g_3D,op.eClosure,id_t)
+            e_l = boundary_restriction(g_3D,op.eClosure,id_l)
+            e_r = boundary_restriction(g_3D,op.eClosure,id_r)
+            e_s = boundary_restriction(g_3D,op.eClosure,id_s)
+            e_n = boundary_restriction(g_3D,op.eClosure,id_n)
+            e_b = boundary_restriction(g_3D,op.eClosure,id_b)
+            e_t = boundary_restriction(g_3D,op.eClosure,id_t)
             e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r),
                           Pair(id_s,e_s),Pair(id_n,e_n),
                           Pair(id_b,e_b),Pair(id_t,e_t))
 
-            d_l = NormalDerivative(g_3D,op.dClosure,id_l)
-            d_r = NormalDerivative(g_3D,op.dClosure,id_r)
-            d_s = NormalDerivative(g_3D,op.dClosure,id_s)
-            d_n = NormalDerivative(g_3D,op.dClosure,id_n)
-            d_b = NormalDerivative(g_3D,op.dClosure,id_b)
-            d_t = NormalDerivative(g_3D,op.dClosure,id_t)
+            d_l = normal_derivative(g_3D,op.dClosure,id_l)
+            d_r = normal_derivative(g_3D,op.dClosure,id_r)
+            d_s = normal_derivative(g_3D,op.dClosure,id_s)
+            d_n = normal_derivative(g_3D,op.dClosure,id_n)
+            d_b = normal_derivative(g_3D,op.dClosure,id_b)
+            d_t = normal_derivative(g_3D,op.dClosure,id_t)
             d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r),
                           Pair(id_s,d_s),Pair(id_n,d_n),
                           Pair(id_b,d_b),Pair(id_t,d_t))
 
-            H_l = boundary_quadrature(g_3D,op.quadratureClosure,id_r)
-            H_r = boundary_quadrature(g_3D,op.quadratureClosure,id_r)
-            H_s = boundary_quadrature(g_3D,op.quadratureClosure,id_s)
-            H_n = boundary_quadrature(g_3D,op.quadratureClosure,id_n)
-            H_b = boundary_quadrature(g_3D,op.quadratureClosure,id_b)
-            H_t = boundary_quadrature(g_3D,op.quadratureClosure,id_t)
+            H_l = inner_product(boundary_grid(g_3D,id_l),op.quadratureClosure)
+            H_r = inner_product(boundary_grid(g_3D,id_r),op.quadratureClosure)
+            H_s = inner_product(boundary_grid(g_3D,id_s),op.quadratureClosure)
+            H_n = inner_product(boundary_grid(g_3D,id_n),op.quadratureClosure)
+            H_b = inner_product(boundary_grid(g_3D,id_b),op.quadratureClosure)
+            H_t = inner_product(boundary_grid(g_3D,id_t),op.quadratureClosure)
             Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r),
                           Pair(id_s,H_s),Pair(id_n,H_n),
                           Pair(id_b,H_b),Pair(id_t,H_t))
@@ -426,14 +429,15 @@
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
             L = laplace(g_1D, op.innerStencil, op.closureStencils)
-            @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils)
+            @test L == second_derivative(g_1D, op.innerStencil, op.closureStencils)
             @test L isa TensorMapping{T,1,1}  where T
         end
         @testset "3D" begin
             L = laplace(g_3D, op.innerStencil, op.closureStencils)
-            Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1)
-            Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2)
-            Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3)
+            @test L isa TensorMapping{T,3,3} where T
+            Dxx = second_derivative(g_3D, op.innerStencil, op.closureStencils,1)
+            Dyy = second_derivative(g_3D, op.innerStencil, op.closureStencils,2)
+            Dzz = second_derivative(g_3D, op.innerStencil, op.closureStencils,3)
             @test L == Dxx + Dyy + Dzz
             @test L isa TensorMapping{T,3,3} where T
         end
@@ -470,7 +474,8 @@
         # 2nd order interior stencil, 1st order boundary stencil,
         # implies that L*v should be exact for binomials up to order 2.
         @testset "2nd order" begin
-            L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            L = laplace(g_3D,op.innerStencil,op.closureStencils)
             @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
             @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
             @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
@@ -480,7 +485,8 @@
         # 4th order interior stencil, 2nd order boundary stencil,
         # implies that L*v should be exact for binomials up to order 3.
         @testset "4th order" begin
-            L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            L = laplace(g_3D,op.innerStencil,op.closureStencils)
             # NOTE: high tolerances for checking the "exact" differentiation
             # due to accumulation of round-off errors/cancellation errors?
             @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
@@ -492,7 +498,7 @@
     end
 end
 
-@testset "Quadrature diagonal" begin
+@testset "Diagonal-stencil inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
     Lz = 1.
@@ -500,65 +506,37 @@
     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)
-    @testset "quadrature" begin
+    @testset "inner_product" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "0D" begin
+            H = inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
+            @test H == IdentityMapping{Float64}()
+            @test H isa TensorMapping{T,0,0} where T
+        end
         @testset "1D" begin
-            H = quadrature(g_1D,op.quadratureClosure)
-            inner_stencil = Stencil((1.,),center=1)
-            @test H == quadrature(g_1D,inner_stencil,op.quadratureClosure)
+            H = inner_product(g_1D,op.quadratureClosure)
+            inner_stencil = CenteredStencil(1.)
+            @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
             @test H isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            H = quadrature(g_2D,op.quadratureClosure)
-            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
-            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
+            H = inner_product(g_2D,op.quadratureClosure)
+            H_x = inner_product(restrict(g_2D,1),op.quadratureClosure)
+            H_y = inner_product(restrict(g_2D,2),op.quadratureClosure)
             @test H == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
     end
 
-    @testset "boundary_quadrature" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        @testset "1D" begin
-            (id_l, id_r) = boundary_identifiers(g_1D)
-            @test boundary_quadrature(g_1D,op.quadratureClosure,id_l) == IdentityMapping{Float64}()
-            @test boundary_quadrature(g_1D,op.quadratureClosure,id_r) == IdentityMapping{Float64}()
-
-        end
-        @testset "2D" begin
-            (id_w, id_e, id_s, id_n) = boundary_identifiers(g_2D)
-            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
-            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_w) == H_y
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_e) == H_y
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_s) == H_x
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_n) == H_x
-        end
-        @testset "3D" begin
-            (id_w, id_e,
-             id_s, id_n,
-             id_t, id_b) = boundary_identifiers(g_3D)
-            H_xy = quadrature(restrict(g_3D,[1,2]),op.quadratureClosure)
-            H_xz = quadrature(restrict(g_3D,[1,3]),op.quadratureClosure)
-            H_yz = quadrature(restrict(g_3D,[2,3]),op.quadratureClosure)
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_w) == H_yz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_e) == H_yz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_s) == H_xz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_n) == H_xz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_t) == H_xy
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_b) == H_xy
-        end
-    end
-
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            H = quadrature(g_1D,op.quadratureClosure)
+            H = inner_product(g_1D,op.quadratureClosure)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = quadrature(g_2D,op.quadratureClosure)
+            H = inner_product(g_2D,op.quadratureClosure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -575,7 +553,7 @@
 
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -584,7 +562,7 @@
 
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -598,13 +576,13 @@
             u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-4
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-8
             end
@@ -612,27 +590,32 @@
     end
 end
 
-@testset "InverseDiagonalQuadrature" 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))
-    @testset "Constructors" begin
+    @testset "inverse_inner_product" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "0D" begin
+            Hi = inverse_inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
+            @test Hi == IdentityMapping{Float64}()
+            @test Hi isa TensorMapping{T,0,0} where T
+        end
         @testset "1D" begin
-            Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure);
-            inner_stencil = Stencil((1.,),center=1)
+            Hi = inverse_inner_product(g_1D, op.quadratureClosure);
+            inner_stencil = CenteredStencil(1.)
             closures = ()
             for i = 1:length(op.quadratureClosure)
                 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
             end
-            @test Hi == InverseQuadrature(g_1D,inner_stencil,closures)
+            @test Hi == inverse_inner_product(g_1D,closures,inner_stencil)
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
-            Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
-            Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+            Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure)
+            Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure)
             @test Hi == Hi_x⊗Hi_y
             @test Hi isa TensorMapping{T,2,2} where T
         end
@@ -641,12 +624,12 @@
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_1D,op.quadratureClosure)
             @test domain_size(Hi) == size(g_1D)
             @test range_size(Hi) == size(g_1D)
         end
         @testset "2D" begin
-            Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -658,15 +641,15 @@
             u = evalOn(g_1D,x->x^3-x^2+1)
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_1D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_1D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -676,15 +659,15 @@
             u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_2D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_2D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -705,7 +688,7 @@
             @test op_l isa TensorMapping{T,0,1} where T
 
             op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
-            @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper())
+            @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper())
             @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
             @test op_r isa TensorMapping{T,0,1} where T
         end
@@ -836,28 +819,28 @@
 
 end
 
-@testset "BoundaryRestriction" begin
+@testset "boundary_restriction" 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 "boundary_restriction" begin
         @testset "1D" begin
-            e_l = BoundaryRestriction(g_1D,op.eClosure,Lower())
-            @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
+            e_l = boundary_restriction(g_1D,op.eClosure,Lower())
+            @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
             @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
             @test e_l isa BoundaryOperator{T,Lower} where T
             @test e_l isa TensorMapping{T,0,1} where T
 
-            e_r = BoundaryRestriction(g_1D,op.eClosure,Upper())
-            @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
+            e_r = boundary_restriction(g_1D,op.eClosure,Upper())
+            @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
             @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper())
             @test e_r isa BoundaryOperator{T,Upper} where T
             @test e_r isa TensorMapping{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            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
@@ -865,8 +848,8 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+            e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
+            e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
 
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
@@ -877,10 +860,10 @@
         end
 
         @testset "2D" begin
-            e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
-            e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
-            e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,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}())
 
             v = rand(11, 15)
             u = fill(3.124)
@@ -893,25 +876,25 @@
     end
 end
 
-@testset "NormalDerivative" begin
+@testset "normal_derivative" begin
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
-    @testset "Constructors" begin
+    @testset "normal_derivative" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            d_l = NormalDerivative(g_1D, op.dClosure, Lower())
-            @test d_l == NormalDerivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_1D, op.dClosure, Lower())
+            @test d_l == normal_derivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
             @test d_l isa BoundaryOperator{T,Lower} where T
             @test d_l isa TensorMapping{T,0,1} where T
         end
         @testset "2D" begin
             op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
             Ix = IdentityMapping{Float64}((size(g_2D)[1],))
             Iy = IdentityMapping{Float64}((size(g_2D)[2],))
-            d_l = NormalDerivative(restrict(g_2D,1),op.dClosure,Lower())
-            d_r = NormalDerivative(restrict(g_2D,2),op.dClosure,Upper())
+            d_l = normal_derivative(restrict(g_2D,1),op.dClosure,Lower())
+            d_r = normal_derivative(restrict(g_2D,2),op.dClosure,Upper())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa TensorMapping{T,1,2} where T
@@ -925,10 +908,10 @@
         # TODO: Test for higher order polynomials?
         @testset "2nd order" begin
             op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
-            d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
-            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
@@ -938,10 +921,10 @@
 
         @testset "4th order" begin
             op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
-            d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
-            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13