--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Mar 21 13:44:29 2021 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Mar 28 15:48:06 2021 +0200
@@ -1,5 +1,5 @@
 """
-    Laplace{T, Dim, ...} <: TensorMapping{T,Dim,Dim}
+    Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim}
     Laplace(grid::AbstractGrid, fn; order)
 
 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
@@ -9,18 +9,19 @@
 Laplace(grid, fn; order) creates the Laplace operator defined on `grid`,
 where the operators are read from TOML. The differential operator is created
 using `laplace(grid::AbstractGrid,...)`.
+
+Note that all properties of Laplace, excluding the Differential operator `D`, are
+abstract types. For performance reasons, they should therefore be
+accessed via the provided getter functions (e.g `inner_product(::Laplace)`).
+
 """
-struct Laplace{T, Dim, Rb, TMdiffop<:TensorMapping{T,Dim,Dim}, # Differential operator
-                           TMipop<:TensorMapping{T,Dim,Dim}, # Inner product operator
-                           TMbop<:TensorMapping{T,Rb,Dim}, # Boundary operator
-                           TMbqop<:TensorMapping{T,Rb,Rb}, # Boundary quadrature
-                           BID<:BoundaryIdentifier} <: TensorMapping{T,Dim,Dim}
+struct Laplace{T, Dim, TMdiffop<:TensorMapping{T,Dim,Dim}} <: TensorMapping{T,Dim,Dim}
     D::TMdiffop # Differential operator
-    H::TMipop # Inner product operator
-    H_inv::TMipop # Inverse inner product operator
-    e::StaticDict{BID,TMbop} # Boundary restriction operators
-    d::StaticDict{BID,TMbop} # Normal derivative operators
-    H_boundary::StaticDict{BID,TMbqop} # Boundary quadrature operators
+    H::TensorMapping # Inner product operator
+    H_inv::TensorMapping # Inverse inner product operator
+    e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators.
+    d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators
+    H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product?
 end
 export Laplace
 
@@ -50,28 +51,80 @@
     return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs))
 end
 
+
+Base.show(io::IO, L::Laplace{T, Dim}) where {T,Dim} =  print(io, "Laplace{$T, $Dim}(", L.D, L.H, L.H_inv, L.e, L.d, L.H_boundary, ")")
+
 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
 
+
+"""
+    inner_product(L::Lapalace)
+
+Returns the inner product operator associated with `L`
+
+"""
 inner_product(L::Laplace) = L.H
 export inner_product
+
+
+"""
+    inverse_inner_product(L::Lapalace)
+
+Returns the inverse of the inner product operator associated with `L`
+
+"""
 inverse_inner_product(L::Laplace) = L.H_inv
 export inverse_inner_product
+
+
+"""
+    boundary_restriction(L::Lapalace,id::BoundaryIdentifier)
+    boundary_restriction(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
+    boundary_restriction(L::Lapalace,ids...)
+
+Returns boundary restriction operator(s) associated with `L` for the boundary(s)
+identified by id(s).
+
+"""
 boundary_restriction(L::Laplace,id::BoundaryIdentifier) = L.e[id]
 boundary_restriction(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N)
 boundary_restriction(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N)
 export boundary_restriction
+
+
+"""
+    normal_derivative(L::Lapalace,id::BoundaryIdentifier)
+    normal_derivative(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
+    normal_derivative(L::Lapalace,ids...)
+
+Returns normal derivative operator(s) associated with `L` for the boundary(s)
+identified by id(s).
+
+"""
 normal_derivative(L::Laplace,id::BoundaryIdentifier) = L.d[id]
 normal_derivative(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
 normal_derivative(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
 export normal_derivative
+
+
 # TODO: boundary_inner_product?
+"""
+    boundary_quadrature(L::Lapalace,id::BoundaryIdentifier)
+    boundary_quadrature(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
+    boundary_quadrature(L::Lapalace,ids...)
+
+Returns boundary quadrature operator(s) associated with `L` for the boundary(s)
+identified by id(s).
+
+"""
 boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id]
 boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
 boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
 export boundary_quadrature
 
+
 """
     laplace(grid, inner_stencil, closure_stencils)