Mercurial > repos > public > sbplib_julia

changeset 875:067a322e4f73 laplace_benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/laplace_opset
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jan 2022 10:55:08 +0100
parents 7e9ebd572deb (current diff) 6a4d36eccf39 (diff)
children 4f3924293894
files src/SbpOperators/volumeops/laplace/laplace.jl
diffstat 48 files changed, 1709 insertions(+), 803 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,3 @@
+syntax: glob
+docs/build/
+docs/build-local/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,54 @@
+JULIA_DEFAULT=julia --startup-file=no
+JULIA?=$(JULIA_DEFAULT)
+
+# Set the default browser
+WHICH_XDG_OPEN=$(shell which xdg-open)
+WHICH_OPEN=$(shell which open)
+BROWSER_DEFAULT  = $(if $(WHICH_XDG_OPEN), xdg-open)
+BROWSER_DEFAULT := $(if $(BROWSER_DEFAULT), $(BROWSER_DEFAULT), open)
+BROWSER?=$(BROWSER_DEFAULT)
+
+help:
+	@echo 'Targets:'
+	@echo '    help        - Show this help.'
+	@echo '    docs        - Generate docs for webserver deployment.'
+	@echo '    localdocs   - Generate docs for local viewing.'
+	@echo '    opendocs    - Open documentation in the browser remaking it if necessary.'
+	@echo ''
+	@echo 'Variables:'
+	@echo '    JULIA       - Controls which command is used to run julia'
+	@echo '                  Default $(JULIA_DEFAULT)'
+	@echo '    BROWSER     - Sets the command for how to open html files'
+	@echo '                  Default: xdg-open if it exists otherwise open'
+	@echo ''
+	@echo 'Variables can be set on the commandline using the -e flag for make, e.g.'
+	@echo '    make localdocs -e JULIA=path/to/julia'
+	@echo 'or as shell environment variables.'
+
+docs: docs/build
+
+localdocs: docs/build-local
+
+opendocs: localdocs
+	$(BROWSER) docs/build-local/index.html
+
+clean:
+	rm -rf docs/build
+	rm -rf docs/build-local
+
+.PHONY: help docs localdocs opendocs clean
+
+SRC_DIRS = src docs/src
+SRC_FILES_AND_DIRS = $(foreach dir,$(SRC_DIRS),$(shell find $(dir)))
+DEP_IGNORE = %/.DS_Store
+DOCS_DEPENDENCIES = docs/make.jl $(filter-out $(DEP_IGNORE),$(SRC_FILES_AND_DIRS))
+docs/build: $(DOCS_DEPENDENCIES)
+	$(JULIA) --project=docs docs/make.jl --build-dir build --prettyurls
+
+docs/build-local: $(DOCS_DEPENDENCIES)
+	$(JULIA) --project=docs docs/make.jl --build-dir build-local
+
+
+.PHONY: temp
+temp:
+	@echo $(SRC_FILES_AND_DIRS)
--- a/Manifest.toml	Thu Jan 20 21:51:53 2022 +0100
+++ b/Manifest.toml	Thu Jan 27 10:55:08 2022 +0100
@@ -1,41 +1,59 @@
 # This file is machine-generated - editing it directly is not advised
 
-[[Adapt]]
+julia_version = "1.7.0"
+manifest_format = "2.0"
+
+[[deps.Adapt]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 version = "3.3.1"
 
-[[Dates]]
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
-[[Libdl]]
+[[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
-[[LinearAlgebra]]
-deps = ["Libdl"]
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[OffsetArrays]]
+[[deps.OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "2bf78c5fd7fa56d2bbf1efbadd45c1b8789e6f57"
+git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.10.2"
+version = "1.10.8"
 
-[[Printf]]
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
-[[TOML]]
+[[deps.TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 
-[[TiledIteration]]
+[[deps.TiledIteration]]
 deps = ["OffsetArrays"]
-git-tree-sha1 = "52c5f816857bfb3291c7d25420b1f4aca0a74d18"
+git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297"
 uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
-version = "0.3.0"
+version = "0.3.1"
 
-[[Unicode]]
+[[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
--- a/Notes.md	Thu Jan 20 21:51:53 2022 +0100
+++ b/Notes.md	Thu Jan 27 10:55:08 2022 +0100
@@ -53,6 +53,23 @@
 
 * Remove order as a table name and put it as a variable.
 
+### Parsing of stencil sets
+At the moment the only parsing that can be done at the top level is conversion
+from the toml file to a dict of strings. This forces the user to dig through
+the dictionary and apply the correct parsing methods for the different parts,
+e.g. `parse_stencil` or `parse_tuple`. While very flexible there is a tight
+coupling between what is written in the file and what code is run to make data
+in the file usable. While this coupling is hard to avoid it should be made
+explicit. This could be done by putting a reference to a parsing function in
+the operator-storage format or somehow specifying the type of each object.
+This mechanism should be extensible without changing the package. Perhaps
+there could be a way to register parsing functions or object types for the
+toml.
+
+If possible the goal should be for the parsing to get all the way to the
+stencils so that a user calls `read_stencil_set` and gets a
+dictionary-structure containing stencils, tuples, scalars and other types
+ready for input to the methods creating the operators.
 
 ## Variable second derivative
 
@@ -321,3 +338,15 @@
 See [this talk](https://www.youtube.com/watch?v=vPsfZUqI4_0) for some simple ideas for defining effecive memory usage and some comparison with peak performance.
 
 
+## Adjoint as a trait on the sbp_operator level?
+
+It would be nice to have a way of refering to adjoints with resepct to the sbp-inner-product.
+If it was possible you could reduce the number of times you have to deal with the inner product matrix.
+
+Since the LazyOperators package is sort of implementing matrix-free matrices there is no concept of inner products there at the moment. It seems to complicate large parts of the package if this was included there.
+
+A different approach would be to include it as a trait for operators so that you can specify what the adjoint for that operator is.
+
+
+## Name of the `VolumeOperator` type for constant stencils
+It seems that the name is too general. The name of the method `volume_operator` makes sense. It should return different types of `TensorMapping` specialized for the grid. A suggetion for a better name is `ConstantStencilVolumeOperator`
--- a/README.md	Thu Jan 20 21:51:53 2022 +0100
+++ b/README.md	Thu Jan 27 10:55:08 2022 +0100
@@ -27,3 +27,23 @@
 Pkg.test(test_args=["*/lazy_tensor_operations_test.jl", "Grids/*"])
 ```
 will run any file named `lazy_tensor_operations_test.jl` and all the files in the `Grids` folder.
+
+
+## Generating and using the documentation
+Generating the documentation can be done using either `make` or through activating the `docs` environment and including the script `docs/make.jl` at the REPL.
+
+Using `make` there are three targets
+```shell
+make docs
+make localdocs
+make opendocs
+make help
+```
+The first variant generates files suitable for webserver deployment, i.e with `prettyurls=true`. The second generates files sutible for local viewing in a web browser, i.e `prettyurls=false`. To view the documentation locally simply open `docs/build/index.html` in your web browser. The documentation can be automatically built and opened using
+```shell
+make opendocs
+```
+
+When including the `docs/make.jl` script `prettyurls` is set to `false` by default.
+
+Including `docs/make.jl` from the REPL may be preferable when repeatadely building the documentation since this avoids compilation latency.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/Manifest.toml	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,145 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.7.1"
+manifest_format = "2.0"
+
+[[deps.ANSIColoredPrinters]]
+git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
+uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
+version = "0.0.1"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.3.3"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.6"
+
+[[deps.Documenter]]
+deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
+git-tree-sha1 = "f425293f7e0acaf9144de6d731772de156676233"
+uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+version = "0.27.10"
+
+[[deps.IOCapture]]
+deps = ["Logging", "Random"]
+git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
+uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
+version = "0.2.2"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.2"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[deps.OffsetArrays]]
+deps = ["Adapt"]
+git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.10.8"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.Parsers]]
+deps = ["Dates"]
+git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "2.1.3"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.Sbplib]]
+deps = ["TOML", "TiledIteration"]
+path = ".."
+uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
+version = "0.1.0"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.TiledIteration]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297"
+uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
+version = "0.3.1"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/Project.toml	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,3 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Sbplib = "5a373a26-915f-4769-bcab-bf03835de17b"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/make.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,42 @@
+using Documenter
+using Sbplib
+
+using Sbplib.DiffOps
+using Sbplib.Grids
+using Sbplib.LazyTensors
+using Sbplib.RegionIndices
+using Sbplib.SbpOperators
+using Sbplib.StaticDicts
+
+sitename = "Sbplib.jl"
+
+if "--prettyurls" ∈ ARGS
+    prettyurls = true
+else
+    prettyurls = false
+end
+
+if "--build-dir" ∈ ARGS
+    i = findlast(==("--build-dir"), ARGS)
+    build = ARGS[i+1]
+else
+    build = "build-local"
+end
+
+pages = [
+    "Home" => "index.md",
+    "operator_file_format.md",
+    "Submodules" => [
+        "submodules/grids.md",
+        "submodules/diff_ops.md",
+        "submodules/lazy_tensors.md",
+        "submodules/region_indices.md",
+        "submodules/sbp_operators.md",
+        "submodules/static_dicts.md",
+    ],
+    "doc_index.md",
+]
+# This ordering is not respected by @contents. See https://github.com/JuliaDocs/Documenter.jl/issues/936
+
+format=Documenter.HTML(;prettyurls)
+makedocs(;sitename, pages, format, build)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/doc_index.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,4 @@
+# Interface index
+
+```@index
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/index.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,9 @@
+# Sbplib.jl
+
+```@contents
+Depth = 1
+```
+
+> Matrix free finite differences methods for all!
+>
+>  -- The mighty Trocatron, ruler of errors and grids
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/operator_file_format.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,72 @@
+# Operator file format
+
+The intention is that Sbplib.jl should be a general and extensible framework
+for working with finite difference methods. It therefore includes a set of
+tools for storing and sharing operator definitions as well as a set of widely
+used operators.
+
+## Using the included operators
+
+Most users will likely access the included operators by simply passing the
+filename of the wanted operator set to the appropriate function.  The location
+of the included stencil sets can be computed using
+[`sbp_operators_path`](@ref).
+```@meta
+# TODO: provide examples of functions to pass the files to
+```
+Advanced user might want to get access to the individual objects of an
+operator file. This can be accomplished using functions such as
+* [`read_stencil_set`](@ref)
+* [`parse_scalar`](@ref)
+* [`parse_stencil`](@ref)
+* [`parse_tuple`](@ref)
+
+When parsing operator objects they are interpreted using `Rational`s and
+possibly have to be converted to a desired type before use. This allows
+preserving maximum accuracy when needed.
+```@meta
+# TBD: "possibly have to be converted to a desired type before use" Is this the case? Can it be fixed?
+```
+
+## File format
+The file format is based on TOML and can be parsed using `TOML.parse`. A file
+can optionally start with a `[meta]` section which can specify things like who
+the author was, a description and how to cite the operators.
+
+After the `[meta]` section one or more stencil sets follow, each one beginning
+with `[[stencil_set]]`. Each stencil set should include descriptors like
+`order`, `name` or `number_of_bondary_points` to make them unique within the
+TOML-file. What descriptors to use are up to the author of the file to decide.
+
+Beyond identifying information the stencil set can contain any valid TOML.
+This data is then parsed by the functions creating specific operators like
+``D_1`` or ``D_2``.
+
+### Numbers
+Number can be represented as regular TOML numbers e.g. `1`, `-0.4` or
+`4.32e-3`. Alternatively they can be represented as strings which allows
+specifying fraction e.g. `"1/2"` or `"0"`.
+
+All numbers are accurately converted to `Rational`s when using the
+[`parse_scalar`](@ref) function.
+
+### Stencils
+Stencils are parsed using [`parse_stencil`](@ref). They can be specified
+either as a simple arrays
+```toml
+stencil = ["-1/2","0", "1/2"]
+```
+which assumes a centered stencil. Or as a TOML inline table
+```toml
+stencil =  {s = ["-24/17", "59/34", "-4/17", "-3/34", "0", "0"], c = 1},
+```
+which allows specifying the center of the stencil using the key `c`.
+
+## Creating your own operator files
+Operator files can be created either to add new variants of existing types of
+operators like ``D_1`` or ``D_2`` or to describe completely new types of
+operators like for example a novel kind of interpolation operator. In the
+second case new parsing functions are also necessary.
+
+The files can then be used to easily test or share different variants of
+operators.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/submodules/diff_ops.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,23 @@
+# DiffOps
+
+## Contents
+```@contents
+Pages = ["diff_ops.md"]
+```
+
+## Index
+```@index
+Pages = ["diff_ops.md"]
+```
+
+## Public interface
+```@autodocs
+Modules = [DiffOps]
+Private = false # Hide unexported objects
+```
+
+## Internal interface
+```@autodocs
+Modules = [DiffOps]
+Public = false # Hide exported objects
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/submodules/grids.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,23 @@
+# Grids
+
+## Contents
+```@contents
+Pages = ["grids.md"]
+```
+
+## Index
+```@index
+Pages = ["grids.md"]
+```
+
+## Public interface
+```@autodocs
+Modules = [Grids]
+Private = false # Hide unexported objects
+```
+
+## Internal interface
+```@autodocs
+Modules = [Grids]
+Public = false # Hide exported objects
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/submodules/lazy_tensors.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,23 @@
+# LazyTensors
+
+## Contents
+```@contents
+Pages = ["lazy_tensors.md"]
+```
+
+## Index
+```@index
+Pages = ["lazy_tensors.md"]
+```
+
+## Public interface
+```@autodocs
+Modules = [LazyTensors]
+Private = false # Hide unexported objects
+```
+
+## Internal interface
+```@autodocs
+Modules = [LazyTensors]
+Public = false # Hide exported objects
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/submodules/region_indices.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,23 @@
+# RegionIndices
+
+## Contents
+```@contents
+Pages = ["region_indices.md"]
+```
+
+## Index
+```@index
+Pages = ["region_indices.md"]
+```
+
+## Public interface
+```@autodocs
+Modules = [RegionIndices]
+Private = false # Hide unexported objects
+```
+
+## Internal interface
+```@autodocs
+Modules = [RegionIndices]
+Public = false # Hide exported objects
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/submodules/sbp_operators.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,23 @@
+# SbpOperators
+
+## Contents
+```@contents
+Pages = ["sbp_operators.md"]
+```
+
+## Index
+```@index
+Pages = ["sbp_operators.md"]
+```
+
+## Public interface
+```@autodocs
+Modules = [SbpOperators]
+Private = false # Hide unexported objects
+```
+
+## Internal interface
+```@autodocs
+Modules = [SbpOperators]
+Public = false # Hide exported objects
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/submodules/static_dicts.md	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,23 @@
+# StaticDicts
+
+## Contents
+```@contents
+Pages = ["static_dicts.md"]
+```
+
+## Index
+```@index
+Pages = ["static_dicts.md"]
+```
+
+## Public interface
+```@autodocs
+Modules = [StaticDicts]
+Private = false # Hide unexported objects
+```
+
+## Internal interface
+```@autodocs
+Modules = [StaticDicts]
+Public = false # Hide exported objects
+```
--- a/src/DiffOps/DiffOps.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/DiffOps/DiffOps.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -93,6 +93,7 @@
 
 
 """
+    BoundaryCondition
 A BoundaryCondition should implement the method
     sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
 """
--- a/src/LazyTensors/lazy_array.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/LazyTensors/lazy_array.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -42,10 +42,10 @@
 
 """
     LazyElementwiseOperation{T,D,Op} <: LazyArray{T,D}
-Struct allowing for lazy evaluation of elementwise operations on AbstractArrays.
+Struct allowing for lazy evaluation of elementwise operations on `AbstractArray`s.
 
-A LazyElementwiseOperation contains two arrays together with an operation.
-The operations are carried out when the LazyElementwiseOperation is indexed.
+A `LazyElementwiseOperation` contains two arrays together with an operation.
+The operations are carried out when the `LazyElementwiseOperation` is indexed.
 """
 struct LazyElementwiseOperation{T,D,Op,T1<:AbstractArray{T,D},T2<:AbstractArray{T,D}} <: LazyArray{T,D}
     a::T1
--- a/src/LazyTensors/lazy_tensor_operations.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -76,7 +76,7 @@
 """
     TensorMappingComposition{T,R,K,D}
 
-Lazily compose two TensorMappings, so that they can be handled as a single TensorMapping.
+Lazily compose two `TensorMapping`s, so that they can be handled as a single `TensorMapping`.
 """
 struct TensorMappingComposition{T,R,K,D, TM1<:TensorMapping{T,R,K}, TM2<:TensorMapping{T,K,D}} <: TensorMapping{T,R,D}
     t1::TM1
@@ -165,8 +165,8 @@
 apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...]
 
 """
-Base.:∘(tm, tmi)
-Base.:∘(tmi, tm)
+    Base.:∘(tm, tmi)
+    Base.:∘(tmi, tm)
 
 Composes a `Tensormapping` `tm` with an `IdentityMapping` `tmi`, by returning `tm`
 """
@@ -316,7 +316,7 @@
     slice_tuple(t, Val(l), Val(u))
 
 Get a slice of a tuple in a type stable way.
-Equivalent to t[l:u] but type stable.
+Equivalent to `t[l:u]` but type stable.
 """
 function slice_tuple(t,::Val{L},::Val{U}) where {L,U}
     return ntuple(i->t[i+L-1], U-L+1)
@@ -327,7 +327,7 @@
 
 Split the tuple `t` into two parts. the first part is `M` long.
 E.g
-```
+```julia
 split_tuple((1,2,3,4),Val(3)) -> (1,2,3), (4,)
 ```
 """
@@ -357,17 +357,19 @@
 flatten_tuple(t::Tuple) = ((flatten_tuple.(t)...)...,) # simplify?
 flatten_tuple(ts::Vararg) = flatten_tuple(ts)
 
-"""
-   LazyOuterProduct(tms...)
+@doc raw"""
+    LazyOuterProduct(tms...)
 
 Creates a `TensorMappingComposition` for the outerproduct of `tms...`.
 This is done by separating the outer product into regular products of outer products involving only identity mappings and one non-identity mapping.
 
 First let
 ```math
-A = A_{I,J}
-B = B_{M,N}
-C = C_{P,Q}
+\begin{aligned}
+A &= A_{I,J} \\
+B &= B_{M,N} \\
+C &= C_{P,Q} \\
+\end{aligned}
 ```
 
 where ``I``, ``M``, ``P`` are  multi-indexes for the ranges of ``A``, ``B``, ``C``, and ``J``, ``N``, ``Q`` are multi-indexes of the domains.
@@ -385,10 +387,12 @@
 ```math
 A⊗B⊗C = (A⊗I_{|M|}⊗I_{|P|})(I_{|J|}⊗B⊗I_{|P|})(I_{|J|}⊗I_{|N|}⊗C)
 ```
-where |.| of a multi-index is a vector of sizes for each dimension. ``I_v`` denotes the identity tensor of size ``v[i]`` in each direction
+where ``|⋅|`` of a multi-index is a vector of sizes for each dimension. ``I_v`` denotes the identity tensor of size ``v[i]`` in each direction
 To apply ``A⊗B⊗C`` we evaluate
 
+```math
 (A⊗B⊗C)v = [(A⊗I_{|M|}⊗I_{|P|})  [(I_{|J|}⊗B⊗I_{|P|}) [(I_{|J|}⊗I_{|N|}⊗C)v]]]
+```
 """
 function LazyOuterProduct end
 export LazyOuterProduct
--- a/src/LazyTensors/tensor_mapping.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/LazyTensors/tensor_mapping.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,21 +1,24 @@
 """
     TensorMapping{T,R,D}
 
-Describes a mapping of a D dimension tensor to an R dimension tensor.
+Describes a mapping of a `D` dimension tensor to an `R` dimension tensor.
 The action of the mapping is implemented through the method
-
+```julia
     apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+```
 
 The size of the range and domain that the operator works with should be returned by
 the functions
-
+```julia
     range_size(::TensorMapping)
     domain_size(::TensorMapping)
-
+```
 to allow querying for one or the other.
 
 Optionally the action of the transpose may be defined through
+```julia
     apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+```
 """
 abstract type TensorMapping{T,R,D} end
 export TensorMapping
@@ -37,11 +40,13 @@
 export apply_transpose
 
 """
+    range_dim(::TensorMapping)
 Return the dimension of the range space of a given mapping
 """
 range_dim(::TensorMapping{T,R,D}) where {T,R,D} = R
 
 """
+    domain_dim(::TensorMapping)
 Return the dimension of the domain space of a given mapping
 """
 domain_dim(::TensorMapping{T,R,D}) where {T,R,D} = D
--- a/src/SbpOperators/SbpOperators.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -5,11 +5,16 @@
 using Sbplib.Grids
 using Sbplib.StaticDicts
 
+@enum Parity begin
+    odd = -1
+    even = 1
+end
+
 include("stencil.jl")
-include("d2.jl")
 include("readoperator.jl")
 include("volumeops/volume_operator.jl")
-include("volumeops/derivatives/secondderivative.jl")
+include("volumeops/constant_interior_scaling_operator.jl")
+include("volumeops/derivatives/second_derivative.jl")
 include("volumeops/laplace/laplace.jl")
 include("volumeops/inner_products/inner_product.jl")
 include("volumeops/inner_products/inverse_inner_product.jl")
@@ -17,4 +22,11 @@
 include("boundaryops/boundary_restriction.jl")
 include("boundaryops/normal_derivative.jl")
 
+
+export inner_product
+export inverse_inner_product
+export boundary_restriction
+export normal_derivative
+export boundary_quadrature
+
 end # module
--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -9,7 +9,7 @@
 of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
 the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
 """
-function boundary_operator(grid::EquidistantGrid{Dim,T}, closure_stencil::Stencil{T}, boundary::CartesianBoundary) where {Dim,T}
+function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
     #TODO:Check that dim(boundary) <= Dim?
 
     # Create 1D boundary operator
@@ -18,8 +18,8 @@
     op = BoundaryOperator(restrict(grid, d), closure_stencil, r)
 
     # Create 1D IdentityMappings for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
-    Is = IdentityMapping{T}.(size.(one_d_grids))
+    one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
+    Is = IdentityMapping{eltype(grid)}.(size.(one_d_grids))
 
     # Formulate the correct outer product sequence of the identity mappings and
     # the boundary operator
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -9,7 +9,8 @@
 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.
 """
-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 boundary_restriction
+function boundary_restriction(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
+    converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil)
+    return SbpOperators.boundary_operator(grid, converted_stencil, boundary)
+end
+boundary_restriction(grid::EquidistantGrid{1}, closure_stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -9,10 +9,9 @@
 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 normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+function normal_derivative(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
 end
-normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
-export normal_derivative
+normal_derivative(grid::EquidistantGrid{1}, closure_stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
--- a/src/SbpOperators/d2.jl	Thu Jan 20 21:51:53 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-export D2, closuresize
-
-@enum Parity begin
-    odd = -1
-    even = 1
-end
-
-struct D2{T,M}
-    innerStencil::Stencil{T}
-    closureStencils::NTuple{M,Stencil{T}}
-    eClosure::Stencil{T}
-    dClosure::Stencil{T}
-    quadratureClosure::NTuple{M,Stencil{T}}
-    parity::Parity
-end
-
-closuresize(D::D2{T,M}) where {T,M} = M
--- a/src/SbpOperators/operators/standard_diagonal.toml	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Thu Jan 27 10:55:08 2022 +0100
@@ -1,36 +1,60 @@
 [meta]
 authors = "Ken Mattson"
-descripion = "Standard operators for equidistant grids"
+description = "Standard operators for equidistant grids"
 type = "equidistant"
+cite = """
+    Ken Mattsson, Jan Nordström,
+    Summation by parts operators for finite difference approximations of second derivatives,
+    Journal of Computational Physics,
+    Volume 199, Issue 2,
+    2004,
+    Pages 503-540,
+    ISSN 0021-9991,
+    https://doi.org/10.1016/j.jcp.2004.03.001.
+"""
 
-[order2]
-H.inner = ["1"]
+[[stencil_set]]
+
+order = 2
+
+H.inner = "1"
 H.closure = ["1/2"]
 
 D1.inner_stencil = ["-1/2", "0", "1/2"]
 D1.closure_stencils = [
-    ["-1", "1"],
+    {s = ["-1", "1"], c = 1},
 ]
 
 D2.inner_stencil = ["1", "-2", "1"]
 D2.closure_stencils = [
-    ["1", "-2", "1"],
+    {s = ["1", "-2", "1"], c = 1},
 ]
 
 e.closure = ["1"]
-d1.closure = ["-3/2", "2", "-1/2"]
+d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1}
+
+[[stencil_set]]
+
+order = 4
 
-[order4]
-H.inner = ["1"]
+H.inner = "1"
 H.closure = ["17/48", "59/48", "43/48", "49/48"]
 
+D1.inner_stencil = ["1/12","-2/3","0","2/3","-1/12"]
+D1.closure_stencils = [
+    {s = [ "-24/17",  "59/34",  "-4/17", "-3/34",     "0",     "0"], c = 1},
+    {s = [   "-1/2",      "0",    "1/2",     "0",     "0",     "0"], c = 2},
+    {s = [   "4/43", "-59/86",      "0", "59/86", "-4/43",     "0"], c = 3},
+    {s = [   "3/98",      "0", "-59/98",     "0", "32/49", "-4/49"], c = 4},
+]
+
 D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
 D2.closure_stencils = [
-    [     "2",    "-5",      "4",       "-1",     "0",     "0"],
-    [     "1",    "-2",      "1",        "0",     "0",     "0"],
-    [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"],
-    [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"],
+    {s = [     "2",    "-5",      "4",       "-1",     "0",     "0"], c = 1},
+    {s = [     "1",    "-2",      "1",        "0",     "0",     "0"], c = 2},
+    {s = [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"], c = 3},
+    {s = [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"], c = 4},
 ]
 
 e.closure = ["1"]
-d1.closure = ["-11/6", "3", "-3/2", "1/3"]
+d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1}
--- a/src/SbpOperators/readoperator.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/readoperator.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,155 +1,160 @@
 using TOML
 
-export read_D2_operator
-export read_stencil
-export read_stencils
-export read_tuple
+export read_stencil_set
+export get_stencil_set
+
+export parse_stencil
+export parse_scalar
+export parse_tuple
+
+export sbp_operators_path
+
 
-export get_stencil
-export get_stencils
-export get_tuple
+"""
+    read_stencil_set(filename; filters)
 
-function read_D2_operator(fn; order)
-    operators = TOML.parsefile(fn)["order$order"]
-    D2 = operators["D2"]
-    H = operators["H"]
-    e = operators["e"]
-    d1 = operators["d1"]
+Picks out a stencil set from a TOML file based on some key-value
+filters. If more than one set matches the filters an error is raised. The
+returned stencil set contains parsed TOML intended for functions like
+`parse_scalar` and `parse_stencil`.
+
+The stencil set is not parsed beyond the inital TOML parse. To get usable
+stencils use the `parse_stencil` functions on the fields of the stencil set.
+
+The reason for this is that since stencil sets are intended to be very
+general, and currently do not include any way to specify how to parse a given
+section, the exact parsing is left to the user.
 
-    # Create inner stencil
-    innerStencil = get_stencil(operators, "D2", "inner_stencil")
+For more information see [Operator file format](@ref) in the documentation.
+
+See also [`sbp_operators_path`](@ref), [`get_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref),.
+"""
+read_stencil_set(filename; filters...) = get_stencil_set(TOML.parsefile(filename); filters...)
+
+"""
+    get_stencil_set(parsed_toml; filters...)
+
+Picks out a stencil set from an already parsed TOML based on some key-value
+filters.
 
-    # Create boundary stencils
-    boundarySize = length(D2["closure_stencils"])
-    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
-    for i ∈ 1:boundarySize
-        closureStencils = (closureStencils..., get_stencil(operators, "D2", "closure_stencils", i; center=i))
+See also [`read_stencil_set`](@ref).
+"""
+function get_stencil_set(parsed_toml; filters...)
+    matches = findall(parsed_toml["stencil_set"]) do set
+        for (key, val) ∈ filters
+            if set[string(key)] != val
+                return false
+            end
+        end
+
+        return true
     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)
 
-    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))
+    if length(matches) != 1
+        throw(ArgumentError("filters must pick out a single stencil set"))
     end
 
-    d2 = SbpOperators.D2(
-        innerStencil,
-        closureStencils,
-        eClosure,
-        dClosure,
-        quadratureClosure,
-        even
-    )
+    i = matches[1]
+    return parsed_toml["stencil_set"][i]
+end
+
+"""
+    parse_stencil(parsed_toml)
+
+Accepts parsed TOML and reads it as a stencil.
+
+See also [`read_stencil_set`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref).
+"""
+function parse_stencil(parsed_toml)
+    check_stencil_toml(parsed_toml)
+
+    if parsed_toml isa Array
+        weights = parse_rational.(parsed_toml)
+        return CenteredStencil(weights...)
+    end
+
+    weights = parse_rational.(parsed_toml["s"])
+    return Stencil(weights..., center = parsed_toml["c"])
+end
+
+"""
+    parse_stencil(T, parsed_toml)
+
+Parses the input as a stencil with element type `T`.
+"""
+parse_stencil(T, parsed_toml) = Stencil{T}(parse_stencil(parsed_toml))
+
+function check_stencil_toml(parsed_toml)
+    if !(parsed_toml isa Dict || parsed_toml isa Vector{String})
+        throw(ArgumentError("the TOML for a stencil must be a vector of strings or a table."))
+    end
+
+    if parsed_toml isa Vector{String}
+        return
+    end
 
-    return d2
+    if !(haskey(parsed_toml, "s") && haskey(parsed_toml, "c"))
+        throw(ArgumentError("the table form of a stencil must have fields `s` and `c`."))
+    end
+
+    if !(parsed_toml["s"] isa Vector{String})
+        throw(ArgumentError("a stencil must be specified as a vector of strings."))
+    end
+
+    if !(parsed_toml["c"] isa Int)
+        throw(ArgumentError("the center of a stencil must be specified as an integer."))
+    end
+end
+
+"""
+    parse_scalar(parsed_toml)
+
+Parse a scalar, represented as a string or a number in the TOML, and return it as a `Rational`
+
+See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref) [`parse_tuple`](@ref).
+"""
+function parse_scalar(parsed_toml)
+    try
+        return parse_rational(parsed_toml)
+    catch e
+        throw(ArgumentError("must be a number or a string representing a number."))
+    end
+end
+
+"""
+    parse_tuple(parsed_toml)
+
+Parse an array as a tuple of scalars.
+
+See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref).
+"""
+function parse_tuple(parsed_toml)
+    if !(parsed_toml isa Array)
+        throw(ArgumentError("argument must be an array"))
+    end
+    return Tuple(parse_scalar.(parsed_toml))
 end
 
 
 """
-    read_stencil(fn, path...; [center])
-
-Read a stencil at `path` from the file with name `fn`.
-If a center is specified the given element of the stecil is set as the center.
-
-See also: [`read_stencils`](@ref), [`read_tuple`](@ref), [`get_stencil`](@ref).
+    parse_rational(parsed_toml)
 
-# Examples
-```
-read_stencil(sbp_operators_path()*"standard_diagonal.toml", "order2", "D2", "inner_stencil")
-read_stencil(sbp_operators_path()*"standard_diagonal.toml", "order2", "d1", "closure"; center=1)
-```
-"""
-read_stencil(fn, path...; center=nothing) = get_stencil(TOML.parsefile(fn), path...; center=center)
-
-"""
-    read_stencils(fn, path...; centers)
-
-Read stencils at `path` from the file `fn`.
-Centers of the stencils are specified as a tuple or array in `centers`.
-
-See also: [`read_stencil`](@ref), [`read_tuple`](@ref), [`get_stencils`](@ref).
-"""
-read_stencils(fn, path...; centers) = get_stencils(TOML.parsefile(fn), path...; centers=centers)
-
+Parse a string or a number as a rational.
 """
-    read_tuple(fn, path...)
-
-Read tuple at `path` from the file `fn`.
-
-See also: [`read_stencil`](@ref), [`read_stencils`](@ref), [`get_tuple`](@ref).
-"""
-read_tuple(fn, path...) = get_tuple(TOML.parsefile(fn), path...)
-
-"""
-    get_stencil(parsed_toml, path...; center=nothing)
-
-Same as [`read_stencil`](@ref)) but takes already parsed toml.
-"""
-get_stencil(parsed_toml, path...; center=nothing) = get_stencil(parsed_toml[path[1]], path[2:end]...; center=center)
-function get_stencil(parsed_toml; center=nothing)
-    @assert parsed_toml isa Vector{String}
-    stencil_weights = Float64.(parse_rational.(parsed_toml))
-
-    width = length(stencil_weights)
-
-    if isnothing(center)
-        center = div(width,2)+1
+function parse_rational(parsed_toml)
+    if parsed_toml isa String
+        expr = Meta.parse(replace(parsed_toml, "/"=>"//"))
+        return eval(:(Rational($expr)))
+    else
+        return Rational(parsed_toml)
     end
-
-    return Stencil(stencil_weights..., center=center)
 end
 
 """
-    get_stencils(parsed_toml, path...; centers)
-
-Same as [`read_stencils`](@ref)) but takes already parsed toml.
-"""
-get_stencils(parsed_toml, path...; centers) = get_stencils(parsed_toml[path[1]], path[2:end]...; centers=centers)
-function get_stencils(parsed_toml; centers)
-    @assert parsed_toml isa Vector{Vector{String}}
-    @assert length(centers) == length(parsed_toml)
+    sbp_operators_path()
 
-    stencils = ()
-    for i ∈ 1:length(parsed_toml)
-        stencil = get_stencil(parsed_toml[i], center = centers[i])
-        stencils = (stencils..., stencil)
-    end
+Calculate the path for the operators folder with included stencil sets.
 
-    return stencils
-end
-
-"""
-    get_tuple(parsed_toml, path...)
-
-Same as [`read_tuple`](@ref)) but takes already parsed toml.
+See also [`read_stencil_set`](@ref)
 """
-get_tuple(parsed_toml, path...) = get_tuple(parsed_toml[path[1]], path[2:end]...)
-function get_tuple(parsed_toml)
-    @assert parsed_toml isa Vector{String}
-    t = Tuple(Float64.(parse_rational.(parsed_toml)))
-    return t
-end
-
-# TODO: Probably should be deleted once we have gotten rid of read_D2_operator()
-function toml_string_array_to_tuple(::Type{T}, arr::AbstractVector{String}) where T
-    return Tuple(T.(parse_rational.(arr)))
-end
-
-function parse_rational(str)
-    expr = Meta.parse(replace(str, "/"=>"//"))
-    return eval(:(Rational($expr)))
-end
-
-function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
-    if N >= n
-        return t
-    else
-        return pad_tuple((t..., zero(T)), n)
-    end
-end
-
 sbp_operators_path() = (@__DIR__) * "/operators/"
-export sbp_operators_path
--- a/src/SbpOperators/stencil.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/stencil.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,10 +1,10 @@
 export CenteredStencil
 
-struct Stencil{T<:Real,N}
+struct Stencil{T,N}
     range::Tuple{Int,Int}
     weights::NTuple{N,T}
 
-    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
+    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N}
         @assert range[2]-range[1]+1 == N
         new{T,N}(range,weights)
     end
@@ -15,13 +15,19 @@
 
 Create a stencil with the given weights with element `center` as the center of the stencil.
 """
-function Stencil(weights::Vararg{Number}; center::Int)
+function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error
     N = length(weights)
     range = (1, N) .- center
 
     return Stencil(range, weights)
 end
 
+function Stencil{T}(s::Stencil) where T
+    return Stencil(s.range, T.(s.weights))
+end
+
+Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil)
+
 function CenteredStencil(weights::Vararg)
     if iseven(length(weights))
         throw(ArgumentError("a centered stencil must have an odd number of weights."))
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,48 @@
+"""
+    ConstantInteriorScalingOperator{T,N} <: TensorMapping{T,1,1}
+
+A one-dimensional operator scaling a vector. The first and last `N` points are
+scaled with individual weights while all interior points are scaled the same.
+"""
+struct ConstantInteriorScalingOperator{T,N} <: TensorMapping{T,1,1}
+    interior_weight::T
+    closure_weights::NTuple{N,T}
+    size::Int
+
+    function ConstantInteriorScalingOperator(interior_weight::T, closure_weights::NTuple{N,T}, size::Int) where {T,N}
+        if size < 2*length(closure_weights)
+            throw(DomainError(size, "size must be larger that two times the closure size."))
+        end
+
+        return new{T,N}(interior_weight, closure_weights, size)
+    end
+end
+
+function ConstantInteriorScalingOperator(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+    return ConstantInteriorScalingOperator(interior_weight, Tuple(closure_weights), size(grid)[1])
+end
+
+closure_size(::ConstantInteriorScalingOperator{T,N}) where {T,N} = N
+
+LazyTensors.range_size(op::ConstantInteriorScalingOperator) = (op.size,)
+LazyTensors.domain_size(op::ConstantInteriorScalingOperator) = (op.size,)
+
+# TBD: @inbounds in apply methods?
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+    return op.closure_weights[Int(i)]*v[Int(i)]
+end
+
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+    return op.interior_weight*v[Int(i)]
+end
+
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T
+    return op.closure_weights[op.size[1]-Int(i)+1]*v[Int(i)]
+end
+
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T
+    r = getregion(i, closure_size(op), op.size[1])
+    return LazyTensors.apply(op, v, Index(i, r))
+end
+
+LazyTensors.apply_transpose(op::ConstantInteriorScalingOperator, v, i) = apply(op, v, i)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,19 @@
+"""
+    second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+
+Creates the second-derivative operator `D2` as a `TensorMapping`
+
+`D2` approximates the second-derivative d²/dξ² on `grid` 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 outer product of the
+one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
+Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
+"""
+function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+    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
+second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
+export second_derivative
--- a/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Thu Jan 20 21:51:53 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,20 +0,0 @@
-"""
-    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`
-
-`D2` approximates the second-derivative d²/dξ² on `grid` 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 outer product of the
-one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
-Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
-"""
-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
-second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
-export second_derivative
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,29 +1,33 @@
 """
-    inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
+    inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
 
-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`.
+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.
+`inner_product` creates `H` on `grid` using the `interior_weight` for the
+interior points and the `closure_weights` for the points close to the
+boundary.
 
-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`,
+On a 1-dimensional grid, `H` is a `ConstantInteriorScalingOperator`. On a
+N-dimensional grid, `H` is the outer product of the 1-dimensional inner
+product operators for each coordinate direction. 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ᵢ
+function inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
+    Hs = ()
+
+    for i ∈ 1:dimension(grid)
+        Hs = (Hs..., inner_product(restrict(grid, i), interior_weight, closure_weights))
     end
+
+    return foldl(⊗, Hs)
+end
+
+function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+    h = spacing(grid)[1]
+
+    H = SbpOperators.ConstantInteriorScalingOperator(grid, h*interior_weight, h.*closure_weights)
     return H
 end
-export inner_product
 
-inner_product(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()
+inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,43 +1,29 @@
 """
-    inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
-    inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})
+    inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
 
-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.
+Constructs the inverse inner product operator `H⁻¹` as a `TensorMapping` using
+the weights of `H`, `interior_weight`, `closure_weights`. `H⁻¹` is inverse of
+the inner product operator `H`.
 
-`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 `ConstantInteriorScalingOperator`. On an
+N-dimensional grid, `H⁻¹` is the outer product of the 1-dimensional inverse
+inner product operators for each coordinate direction. On a 0-dimensional
+`grid`, `H⁻¹` is a 0-dimensional `IdentityMapping`.
+"""
+function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
+    H⁻¹s = ()
 
-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ᵢ⁻¹
+    for i ∈ 1:dimension(grid)
+        H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), interior_weight, closure_weights))
     end
+
+    return foldl(⊗, H⁻¹s)
+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⁻¹
 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)
+inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,19 +1,39 @@
+export Laplace
+export laplace
+# REVIEW: Makes more sense to me to have the exports at the top of the file.
+# Might as well start fixing that.
+
+# REVIEW: The style of name `Laplace` might clash with other concepts. When
+# thinking about implementing the variable second derivative I think I will
+# have to create it as a full TM for the full dimensional problem instead of
+# building it as a 1D operator and then use that with outer products. The
+# natural name there would be `VariableSecondDerivative` (or something
+# similar). But the similarity of the two names would suggest that `Laplace`
+# and `VariableSecondDerivative` are the same kind of thing, which they
+# wouldn't be.
+#
+# How do we distinguish the kind of type we are implementing here and what we
+# could potentially do for the variable second derivative?
+#
+# I see two ways out:
+#   * Come up with a name for these sets of operators and change `Laplace` accordingly.
+#   * Come up with a name for the bare operators and change `VariableSecondDerivative` accordingly.
+
 """
     Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim}
-    Laplace(grid::AbstractGrid, fn; order)
+    Laplace(grid, filename; order)
 
 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
 `TensorMapping`. Additionally, `Laplace` stores the inner product and boundary
-operators relevant for constructing a SBP finite difference scheme as `TensorMapping`s.
+operators relevant for constructing a SBP finite difference scheme as a `TensorMapping`.
 
-Laplace(grid, fn; order) creates the Laplace operator defined on `grid`,
+`Laplace(grid, filename; order)` creates the Laplace operator defined on `grid`,
 where the operators are read from TOML. The differential operator is created
-using `laplace(grid::AbstractGrid,...)`.
+using `laplace(grid,...)`.
 
-Note that all properties of Laplace, excluding the Differential operator `D`, are
+Note that all properties of Laplace, excluding the differential operator `Laplace.D`, are
 abstract types. For performance reasons, they should therefore be
 accessed via the provided getter functions (e.g `inner_product(::Laplace)`).
-
 """
 struct Laplace{T, Dim, TMdiffop<:TensorMapping{T,Dim,Dim}} <: TensorMapping{T,Dim,Dim}
     D::TMdiffop # Differential operator
@@ -21,39 +41,41 @@
     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?
-    order::Int
+    H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators
 end
-export Laplace
 
-function Laplace(grid::AbstractGrid, fn; order)
+function Laplace(grid, filename; order)
+    
+    # Read stencils
+    stencil_set = read_stencil_set(filename; order)
     # TODO: Removed once we can construct the volume and
-    # boundary operators by op(grid, fn; order,...).
-    # Read stencils
-    op = read_D2_operator(fn; order)
-    D_inner_stecil = op.innerStencil
-    D_closure_stencils = op.closureStencils
-    H_closure_stencils = op.quadratureClosure
-    e_closure_stencil = op.eClosure
-    d_closure_stencil = op.dClosure
+    # boundary operators by op(grid, read_stencil_set(fn; order,...)).
+    D_inner_stecil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+    D_closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+    H_inner_stencils = parse_scalar(stencil_set["H"]["inner"])
+    H_closure_stencils = parse_tuple(stencil_set["H"]["closure"])
+    e_closure_stencil = parse_stencil(stencil_set["e"]["closure"])
+    d_closure_stencil = parse_stencil(stencil_set["d1"]["closure"])
+    # REVIEW: Do we add the methods to get rid of this in this branch or a new one?
 
     # Volume operators
     Δ =  laplace(grid, D_inner_stecil, D_closure_stencils)
-    H =  inner_product(grid, H_closure_stencils)
-    H⁻¹ =  inverse_inner_product(grid, H_closure_stencils)
+    H =  inner_product(grid, H_inner_stencils, H_closure_stencils)
+    H⁻¹ = inverse_inner_product(grid, H_inner_stencils, H_closure_stencils)
 
     # Boundary operator - id pairs
     ids = boundary_identifiers(grid)
-    n_ids = length(ids)
-    e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid,e_closure_stencil,ids[i]),n_ids)
-    d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids)
-    Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids)
+    # REVIEW: Change suggestion: Seems more readable to me but pretty subjective so feel free to ignore
+    e_pairs  = map(id -> Pair(id, boundary_restriction(grid, e_closure_stencil, id)),                            ids)
+    d_pairs  = map(id -> Pair(id, normal_derivative(grid, d_closure_stencil, id)),                               ids)
+    Hᵧ_pairs = map(id -> Pair(id, inner_product(boundary_grid(grid, id), H_inner_stencils, H_closure_stencils)), ids)
 
     return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs), order)
 end
 
 # TODO: Consider pretty printing of the following form
 # Base.show(io::IO, L::Laplace{T, Dim}) where {T,Dim,TM} =  print(io, "Laplace{$T, $Dim, $TM}(", L.D, L.H, L.H_inv, L.e, L.d, L.H_boundary, ")")
+# REVIEW: Should leave a todo here to update this once we have some pretty printing for tensor mappings in general.
 
 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
@@ -69,69 +91,65 @@
 export closure_size
 
 """
-    inner_product(L::Lapalace)
+    inner_product(L::Laplace)
 
 Returns the inner product operator associated with `L`
-
 """
 inner_product(L::Laplace) = L.H
-export inner_product
 
 
 """
-    inverse_inner_product(L::Lapalace)
+    inverse_inner_product(L::Laplace)
 
 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...)
+    boundary_restriction(L::Laplace, id::BoundaryIdentifier)
+    boundary_restriction(L::Laplace, ids::Tuple)
+    boundary_restriction(L::Laplace, ids...)
 
 Returns boundary restriction operator(s) associated with `L` for the boundary(s)
 identified by id(s).
+"""
+boundary_restriction(L::Laplace, id::BoundaryIdentifier) = L.e[id]
+boundary_restriction(L::Laplace, ids::Tuple) = map(id-> L.e[id], ids)
+boundary_restriction(L::Laplace, ids...) = boundary_restriction(L, ids)
+# REVIEW: I propose changing the following implementations according to the
+# above. There are some things we're missing with regards to
+# `BoundaryIdentifier`, for example we should be able to handle groups of
+# boundaries as a single `BoundaryIdentifier`. I don't know if we can figure
+# out the interface for that now or if we save it for the future but either
+# way these methods will be affected.
 
-"""
-boundary_restriction(L::Laplace,id::BoundaryIdentifier) = L.e[id]
-boundary_restriction(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N)
-boundary_restriction(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N)
-export boundary_restriction
 
 
 """
-    normal_derivative(L::Lapalace,id::BoundaryIdentifier)
-    normal_derivative(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
-    normal_derivative(L::Lapalace,ids...)
+    normal_derivative(L::Laplace, id::BoundaryIdentifier)
+    normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
+    normal_derivative(L::Laplace, ids...)
 
 Returns normal derivative operator(s) associated with `L` for the boundary(s)
 identified by id(s).
+"""
+normal_derivative(L::Laplace, id::BoundaryIdentifier) = L.d[id]
+normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
+normal_derivative(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
+
 
 """
-normal_derivative(L::Laplace,id::BoundaryIdentifier) = L.d[id]
-normal_derivative(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
-normal_derivative(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
-export normal_derivative
-
-
-# TODO: boundary_inner_product?
-"""
-    boundary_quadrature(L::Lapalace,id::BoundaryIdentifier)
-    boundary_quadrature(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
-    boundary_quadrature(L::Lapalace,ids...)
+    boundary_quadrature(L::Laplace, id::BoundaryIdentifier)
+    boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
+    boundary_quadrature(L::Laplace, ids...)
 
 Returns boundary quadrature operator(s) associated with `L` for the boundary(s)
 identified by id(s).
-
 """
-boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id]
-boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
-boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
-export boundary_quadrature
+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)
 
 abstract type BoundaryConditionType end
 struct NeumannBC <: BoundaryConditionType end
@@ -153,53 +171,23 @@
     return closure
 end
 
-function dirichlet_bc(L::Laplace, id::BoundaryIdentifier)
-    error("Not implemented")
-    # H_inv = inverse_inner_product(L)
-    # e = boundary_restriction(L,id)
-    # d = normal_derivative(L,id)
-    # H_b = boundary_quadrature(L,id)
-    # gamma = borrowing_parameter(L)
-    # tuning = 1.2
-    # S = ScalingOperator(tuning * -1.0/gamma)
-    # tau = H_inv∘(S∘e' + d')∘H_b
-    # closure = tau∘e
-    # penalty = ScalingOperator(-1)∘tau
-    # return (closure, penalty)
-end
-
-# function borrowing_parameter(L)
-#     if L.order == 2
-#         return 0.4
-#     elseif L.order == 4
-#         return 0.2508
-#     elseif L.order == 6
-#         return 0.1878
-#     elseif L.order == 8
-#         return 0.0015
-#     elseif L.order == 10
-#         return 0.0351
-#     end
-# end
-
 """
-    laplace(grid, inner_stencil, closure_stencils)
+    laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
 
 Creates the Laplace operator operator `Δ` as a `TensorMapping`
 
-`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,N on the N-dimensional
-`grid`, using the stencil `inner_stencil` in the interior and a set of stencils
-`closure_stencils` for the points in the closure regions.
+`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using
+the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
+for the points in the closure regions.
 
 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
 where the sum is carried out lazily.
 """
-function laplace(grid::AbstractGrid, inner_stencil, closure_stencils)
+function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
     Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
     for d = 2:dimension(grid)
         Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
     end
     return Δ
 end
-export laplace
--- a/src/SbpOperators/volumeops/volume_operator.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -1,21 +1,22 @@
 """
-    volume_operator(grid,inner_stencil,closure_stencils,parity,direction)
+    volume_operator(grid, inner_stencil, closure_stencils, parity, direction)
+
 Creates a volume operator on a `Dim`-dimensional grid acting along the
-specified coordinate `direction`. The action of the operator is determined by the
-stencils `inner_stencil` and `closure_stencils`.
-When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned.
-When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product
-of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
-the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
+specified coordinate `direction`. The action of the operator is determined by
+the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the
+corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the
+returned operator is the appropriate outer product of a one-dimensional
+operators and `IdentityMapping`s, e.g for `Dim=3` the volume operator in the
+y-direction is `I⊗op⊗I`.
 """
-function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M}
+function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction)
     #TODO: Check that direction <= Dim?
 
     # Create 1D volume operator in along coordinate direction
     op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity)
     # Create 1D IdentityMappings for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
-    Is = IdentityMapping{T}.(size.(one_d_grids))
+    one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
+    Is = IdentityMapping{eltype(grid)}.(size.(one_d_grids))
     # Formulate the correct outer product sequence of the identity mappings and
     # the volume operator
     parts = Base.setindex(Is, op, direction)
@@ -34,7 +35,7 @@
 end
 
 function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity)
-    return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity)
+    return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity)
 end
 
 closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M
--- a/src/StaticDicts/StaticDicts.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/src/StaticDicts/StaticDicts.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -10,9 +10,9 @@
 
 The immutable nature means that `StaticDict` can be compared with `===`, in
 constrast to regular `Dict` or `ImmutableDict` which can not. (See
-https://github.com/JuliaLang/julia/issues/4648 for details) One important
+<https://github.com/JuliaLang/julia/issues/4648> for details.) One important
 aspect of this is that `StaticDict` can be used in a struct while still
-allowing the struct to be comared using the default implementation of `==` for
+allowing the struct to be compared using the default implementation of `==` for
 structs.
 
 Lookups are done by linear search.
--- a/test/Manifest.toml	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/Manifest.toml	Thu Jan 27 10:55:08 2022 +0100
@@ -1,144 +1,264 @@
 # This file is machine-generated - editing it directly is not advised
 
-[[Artifacts]]
-deps = ["Pkg"]
-git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744"
+julia_version = "1.7.0"
+manifest_format = "2.0"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
-version = "1.3.0"
 
-[[Base64]]
+[[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
-[[CompilerSupportLibraries_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "8e695f735fca77e9708e795eda62afdb869cbb70"
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.11.2"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.2"
+
+[[deps.Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.41.0"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "0.3.4+0"
 
-[[Dates]]
+[[deps.Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
-[[DeepDiffs]]
+[[deps.DeepDiffs]]
 git-tree-sha1 = "9824894295b62a6a4ab6adf1c7bf337b3a9ca34c"
 uuid = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 version = "1.2.0"
 
-[[DiffRules]]
-deps = ["NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "eb0c34204c8410888844ada5359ac8b96292cfd1"
+[[deps.DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[deps.DiffRules]]
+deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "9bc5dac3c8b6706b58ad5ce24cffd9861f07c94f"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.0.1"
+version = "1.9.0"
 
-[[Distributed]]
+[[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
-[[Glob]]
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.6"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[deps.Glob]]
 git-tree-sha1 = "4df9f7e06108728ebf00a0a11edee4b29a482bb2"
 uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
 version = "1.3.0"
 
-[[InteractiveUtils]]
+[[deps.InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
-[[JLLWrappers]]
-git-tree-sha1 = "7cec881362e5b4e367ff0279dd99a06526d51a55"
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.2"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.1"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.1.2"
+version = "1.3.0"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
 
-[[LibGit2]]
-deps = ["Printf"]
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
-[[Libdl]]
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
-[[LinearAlgebra]]
-deps = ["Libdl"]
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[Logging]]
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.6"
+
+[[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
-[[Markdown]]
+[[deps.Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
-[[NaNMath]]
-git-tree-sha1 = "c84c576296d0e2fbb3fc134d3e09086b3ea617cd"
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[deps.NaNMath]]
+git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "0.3.4"
+version = "0.3.6"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 
-[[OpenSpecFun_jll]]
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+
+[[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3"
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.3+4"
+version = "0.5.5+0"
 
-[[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
-[[Printf]]
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.2.2"
+
+[[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
-[[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
-[[Random]]
-deps = ["Serialization"]
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[[Requires]]
+[[deps.Requires]]
 deps = ["UUIDs"]
-git-tree-sha1 = "28faf1c963ca1dc3ec87f166d92982e3c4a1f66d"
+git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.1.0"
+version = "1.2.0"
 
-[[SHA]]
+[[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
-[[Serialization]]
+[[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
-[[Sockets]]
+[[deps.SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
-[[SpecialFunctions]]
-deps = ["OpenSpecFun_jll"]
-git-tree-sha1 = "d8d8b8a9f4119829410ecd706da4cc8594a1e020"
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "0.10.3"
+version = "2.0.0"
 
-[[TOML]]
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.TOML]]
 deps = ["Dates"]
-git-tree-sha1 = "d0ac7eaad0fb9f6ba023a1d743edca974ae637c4"
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
-version = "1.0.0"
 
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[[TestSetExtensions]]
+[[deps.TestSetExtensions]]
 deps = ["DeepDiffs", "Distributed", "Test"]
 git-tree-sha1 = "3a2919a78b04c29a1a57b05e1618e473162b15d0"
 uuid = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 version = "2.0.0"
 
-[[Tullio]]
-deps = ["DiffRules", "LinearAlgebra", "Requires"]
-git-tree-sha1 = "b27ec3ce782f69c1c24f373bfb6aa60300ed57c7"
+[[deps.Tullio]]
+deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"]
+git-tree-sha1 = "0288b7a395fc412952baf756fac94e4f28bfec65"
 uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
-version = "0.2.8"
+version = "0.3.2"
 
-[[UUIDs]]
+[[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
-[[Unicode]]
+[[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -8,27 +8,28 @@
 import Sbplib.SbpOperators.BoundaryOperator
 
 @testset "boundary_restriction" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
+	e_closure = parse_stencil(stencil_set["e"]["closure"])
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            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())
+            e_l = boundary_restriction(g_1D,e_closure,Lower())
+            @test e_l == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
+            @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Lower())
             @test e_l isa BoundaryOperator{T,Lower} where T
             @test e_l isa TensorMapping{T,0,1} where T
 
-            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())
+            e_r = boundary_restriction(g_1D,e_closure,Upper())
+            @test e_r == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
+            @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),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 = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            e_w = boundary_restriction(g_2D,e_closure,CartesianBoundary{1,Upper}())
             @test e_w isa InflatedTensorMapping
             @test e_w isa TensorMapping{T,1,2} where T
         end
@@ -36,8 +37,8 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+            e_l = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Lower}())
+            e_r = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Upper}())
 
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
@@ -48,10 +49,10 @@
         end
 
         @testset "2D" begin
-            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}())
+            e_w = boundary_restriction(g_2D, e_closure, CartesianBoundary{1,Lower}())
+            e_e = boundary_restriction(g_2D, e_closure, CartesianBoundary{1,Upper}())
+            e_s = boundary_restriction(g_2D, e_closure, CartesianBoundary{2,Lower}())
+            e_n = boundary_restriction(g_2D, e_closure, CartesianBoundary{2,Upper}())
 
             v = rand(11, 15)
             u = fill(3.124)
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -11,21 +11,21 @@
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
     @testset "normal_derivative" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    	d_closure = parse_stencil(stencil_set["d1"]["closure"])
         @testset "1D" begin
-            d_l = normal_derivative(g_1D, op.dClosure, Lower())
-            @test d_l == normal_derivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_1D, d_closure, Lower())
+            @test d_l == normal_derivative(g_1D, d_closure, 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 = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
+            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
             Ix = IdentityMapping{Float64}((size(g_2D)[1],))
             Iy = IdentityMapping{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(restrict(g_2D,1),op.dClosure,Lower())
-            d_r = normal_derivative(restrict(g_2D,2),op.dClosure,Upper())
+            d_l = normal_derivative(restrict(g_2D,1),d_closure,Lower())
+            d_r = normal_derivative(restrict(g_2D,2),d_closure,Upper())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa TensorMapping{T,1,2} where T
@@ -38,11 +38,12 @@
         v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x)
         # TODO: Test for higher order polynomials?
         @testset "2nd order" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-            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}())
+        	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+        	d_closure = parse_stencil(stencil_set["d1"]["closure"])
+            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
@@ -51,11 +52,12 @@
         end
 
         @testset "4th order" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            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}())
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+        	d_closure = parse_stencil(stencil_set["d1"]["closure"])
+            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
--- a/test/SbpOperators/readoperator_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/readoperator_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -5,6 +5,153 @@
 
 import Sbplib.SbpOperators.Stencil
 
+@testset "readoperator" begin
+    toml_str = """
+        [meta]
+        authors = "Ken Mattson"
+        description = "Standard operators for equidistant grids"
+        type = "equidistant"
+        cite = "A paper a long time ago in a galaxy far far away."
+
+        [[stencil_set]]
+
+        order = 2
+        test = 2
+
+        H.inner = ["1"]
+        H.closure = ["1/2"]
+
+        D1.inner_stencil = ["-1/2", "0", "1/2"]
+        D1.closure_stencils = [
+            {s = ["-1", "1"], c = 1},
+        ]
+
+        D2.inner_stencil = ["1", "-2", "1"]
+        D2.closure_stencils = [
+            {s = ["1", "-2", "1"], c = 1},
+        ]
+
+        e.closure = ["1"]
+        d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1}
+
+        [[stencil_set]]
+
+        order = 4
+        test = 1
+        H.inner = ["1"]
+        H.closure = ["17/48", "59/48", "43/48", "49/48"]
+
+        D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
+        D2.closure_stencils = [
+            {s = [     "2",    "-5",      "4",       "-1",     "0",     "0"], c = 1},
+            {s = [     "1",    "-2",      "1",        "0",     "0",     "0"], c = 2},
+            {s = [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"], c = 3},
+            {s = [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"], c = 4},
+        ]
+
+        e.closure = ["1"]
+        d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1}
+
+        [[stencil_set]]
+        order = 4
+        test = 2
+
+        H.closure = ["-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"]
+    """
+
+    parsed_toml = TOML.parse(toml_str)
+
+    @testset "get_stencil_set" begin
+        @test get_stencil_set(parsed_toml; order = 2) isa Dict
+        @test get_stencil_set(parsed_toml; order = 2) == parsed_toml["stencil_set"][1]
+        @test get_stencil_set(parsed_toml; test = 1) == parsed_toml["stencil_set"][2]
+        @test get_stencil_set(parsed_toml; order = 4, test = 2) == parsed_toml["stencil_set"][3]
+
+        @test_throws ArgumentError get_stencil_set(parsed_toml; test = 2)
+        @test_throws ArgumentError get_stencil_set(parsed_toml; order = 4)
+    end
+
+    @testset "parse_stencil" begin
+        toml = """
+            s1 = ["-1/12","4/3","-5/2","4/3","-1/12"]
+            s2 = {s = ["2", "-5", "4", "-1", "0", "0"], c = 1}
+            s3 = {s = ["1", "-2", "1", "0", "0", "0"], c = 2}
+            s4 = "not a stencil"
+            s5 = [-1, 4, 3]
+            s6 = {k = ["1", "-2", "1", "0", "0", "0"], c = 2}
+            s7 = {s = [-1, 4, 3], c = 2}
+            s8 = {s = ["1", "-2", "1", "0", "0", "0"], c = [2,2]}
+        """
+
+        @test parse_stencil(TOML.parse(toml)["s1"]) == CenteredStencil(-1//12, 4//3, -5//2, 4//3, -1//12)
+        @test parse_stencil(TOML.parse(toml)["s2"]) == Stencil(2//1, -5//1, 4//1, -1//1, 0//1, 0//1; center=1)
+        @test parse_stencil(TOML.parse(toml)["s3"]) == Stencil(1//1, -2//1, 1//1, 0//1, 0//1, 0//1; center=2)
+
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s4"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s5"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s6"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s7"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s8"])
+
+        stencil_set = get_stencil_set(parsed_toml; order = 4, test = 1)
+
+        @test parse_stencil.(stencil_set["D2"]["closure_stencils"]) == [
+            Stencil(  2//1,  -5//1,     4//1,    -1//1,   0//1,   0//1; center=1),
+            Stencil(  1//1,  -2//1,     1//1,     0//1,   0//1,   0//1; center=2),
+            Stencil(-4//43, 59//43, -110//43,   59//43, -4//43,   0//1; center=3),
+            Stencil(-1//49,   0//1,   59//49, -118//49, 64//49, -4//49; center=4),
+        ]
+
+
+        @test parse_stencil(Float64, TOML.parse(toml)["s1"]) == CenteredStencil(-1/12, 4/3, -5/2, 4/3, -1/12)
+        @test parse_stencil(Float64, TOML.parse(toml)["s2"]) == Stencil(2/1, -5/1, 4/1, -1/1, 0/1, 0/1; center=1)
+        @test parse_stencil(Float64, TOML.parse(toml)["s3"]) == Stencil(1/1, -2/1, 1/1, 0/1, 0/1, 0/1; center=2)
+    end
+
+    @testset "parse_scalar" begin
+        toml = TOML.parse("""
+            a1 = 1
+            a2 = 1.5
+            a3 = 1.0
+            a4 = 10
+            a5 = "1/2"
+            a6 = "1.5"
+
+            e1 = [1,2,3]
+            e2 = "a string value"
+        """)
+
+        @test parse_scalar(toml["a1"]) == 1//1
+        @test parse_scalar(toml["a2"]) == 3//2
+        @test parse_scalar(toml["a3"]) == 1//1
+        @test parse_scalar(toml["a4"]) == 10//1
+        @test parse_scalar(toml["a5"]) == 1//2
+        @test parse_scalar(toml["a6"]) == 3//2
+
+        @test_throws ArgumentError parse_scalar(toml["e1"])
+        @test_throws ArgumentError parse_scalar(toml["e2"])
+    end
+
+    @testset "parse_tuple" begin
+        toml = TOML.parse("""
+            t1 = [1,3,4]
+            t2 = ["1/2","3/4","2/1"]
+
+            e1 = "not a tuple"
+            e2.a="1"
+            e3 = 1
+            e4 = ["1/2","3/4","not a number"]
+        """)
+
+        @test parse_tuple(toml["t1"]) == (1//1,3//1,4//1)
+        @test parse_tuple(toml["t2"]) == (1//2,3//4,2//1)
+
+        @test_throws ArgumentError parse_tuple(toml["e1"])
+        @test_throws ArgumentError parse_tuple(toml["e2"])
+        @test_throws ArgumentError parse_tuple(toml["e3"])
+        @test_throws ArgumentError parse_tuple(toml["e4"])
+    end
+end
 
 @testset "parse_rational" begin
     @test SbpOperators.parse_rational("1") isa Rational
@@ -13,81 +160,13 @@
     @test SbpOperators.parse_rational("1/2") == 1//2
     @test SbpOperators.parse_rational("37/13") isa Rational
     @test SbpOperators.parse_rational("37/13") == 37//13
-end
 
-@testset "readoperator" begin
-    toml_str = """
-        [meta]
-        type = "equidistant"
-
-        [order2]
-        H.inner = ["1"]
-
-        D1.inner_stencil = ["-1/2", "0", "1/2"]
-        D1.closure_stencils = [
-            ["-1", "1"],
-        ]
-
-        d1.closure = ["-3/2", "2", "-1/2"]
-
-        [order4]
-        H.closure = ["17/48", "59/48", "43/48", "49/48"]
-
-        D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
-        D2.closure_stencils = [
-            [     "2",    "-5",      "4",       "-1",     "0",     "0"],
-            [     "1",    "-2",      "1",        "0",     "0",     "0"],
-            [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"],
-            [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"],
-        ]
-    """
-
-    parsed_toml = TOML.parse(toml_str)
-    @testset "get_stencil" begin
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil(-1/2, 0., 1/2, center=2)
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil(-1/2, 0., 1/2; center=1)
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil(-1/2, 0., 1/2; center=3)
-
-        @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil(1.; center=1)
+    @test SbpOperators.parse_rational(0.5) isa Rational
+    @test SbpOperators.parse_rational(0.5) == 1//2
 
-        @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
-        @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
-    end
-
-    @testset "get_stencils" begin
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil(-1., 1., center=1),)
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil(-1., 1., center=2),)
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil(-1., 1., center=2),)
-
-        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
-            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=1),
-            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=1),
-            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=1),
-            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=1),
-        )
+    @test SbpOperators.parse_rational("0.5") isa Rational
+    @test SbpOperators.parse_rational("0.5") == 1//2
 
-        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
-            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=4),
-            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=2),
-            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=3),
-            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=1),
-        )
-
-        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
-            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=1),
-            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=2),
-            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=3),
-            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=4),
-        )
-
-        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
-        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4))
-        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2))
-    end
-
-    @testset "get_tuple" begin
-        @test get_tuple(parsed_toml, "order2", "d1", "closure") == (-3/2, 2, -1/2)
-
-        @test_throws AssertionError get_tuple(parsed_toml, "meta", "type")
-    end
+    @test SbpOperators.parse_rational(2) isa Rational
+    @test SbpOperators.parse_rational(2) == 2//1
 end
--- a/test/SbpOperators/stencil_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/stencil_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -15,4 +15,17 @@
 
     @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
     @test_throws ArgumentError CenteredStencil(1,2,3,4)
+
+    # Changing the type of the weights
+    @test Stencil{Float64}(Stencil(1,2,3,4,5; center=2)) == Stencil(1.,2.,3.,4.,5.; center=2)
+    @test Stencil{Float64}(CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
+    @test Stencil{Int}(Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
+    @test Stencil{Rational}(Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+
+    @testset "convert" begin
+        @test convert(Stencil{Float64}, Stencil(1,2,3,4,5; center=2)) == Stencil(1.,2.,3.,4.,5.; center=2)
+        @test convert(Stencil{Float64}, CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
+        @test convert(Stencil{Int}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
+        @test convert(Stencil{Rational}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+    end
 end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,36 @@
+using Test
+
+using Sbplib.LazyTensors
+using Sbplib.SbpOperators
+import Sbplib.SbpOperators: ConstantInteriorScalingOperator
+using Sbplib.Grids
+
+@testset "ConstantInteriorScalingOperator" begin
+    @test ConstantInteriorScalingOperator(1, (2,3), 10) isa ConstantInteriorScalingOperator{Int,2}
+    @test ConstantInteriorScalingOperator(1., (2.,3.), 10) isa ConstantInteriorScalingOperator{Float64,2}
+
+    a = ConstantInteriorScalingOperator(4, (2,3), 10)
+    v = ones(Int, 10)
+    @test a*v == [2,3,4,4,4,4,4,4,3,2]
+    @test a'*v == [2,3,4,4,4,4,4,4,3,2]
+
+    @test range_size(a) == (10,)
+    @test domain_size(a) == (10,)
+
+
+    a = ConstantInteriorScalingOperator(.5, (.1,.2), 7)
+    v = ones(7)
+
+    @test a*v == [.1,.2,.5,.5,.5,.2,.1]
+    @test a'*v == [.1,.2,.5,.5,.5,.2,.1]
+
+    @test range_size(a) == (7,)
+    @test domain_size(a) == (7,)
+
+    @test_throws DomainError ConstantInteriorScalingOperator(4,(2,3), 3)
+
+    @testset "Grid constructor" begin
+        g = EquidistantGrid(11, 0., 2.)
+        @test ConstantInteriorScalingOperator(g, 3., (.1,.2)) isa ConstantInteriorScalingOperator{Float64}
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,118 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+import Sbplib.SbpOperators.VolumeOperator
+
+@testset "SecondDerivative" begin
+    stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; 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))
+
+    @testset "Constructors" begin
+        @testset "1D" begin
+            Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+            @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils,1)
+            @test Dₓₓ isa VolumeOperator
+        end
+        @testset "2D" begin
+            Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1)
+            D2 = second_derivative(g_1D,inner_stencil,closure_stencils)
+            I = IdentityMapping{Float64}(size(g_2D)[2])
+            @test Dₓₓ == D2⊗I
+            @test Dₓₓ isa TensorMapping{T,2,2} where T
+        end
+    end
+
+    # Exact differentiation is measured point-wise. In other cases
+    # the error is measured in the l2-norm.
+    @testset "Accuracy" begin
+        @testset "1D" begin
+            l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
+            monomials = ()
+            maxOrder = 4;
+            for i = 0:maxOrder-1
+                f_i(x) = 1/factorial(i)*x^i
+                monomials = (monomials...,evalOn(g_1D,f_i))
+            end
+            v = evalOn(g_1D,x -> sin(x))
+            vₓₓ = evalOn(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.
+            @testset "2nd order" begin
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+			    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+                Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+                @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
+                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
+            end
+
+            # 4th order interior stencil, 2nd order boundary stencil,
+            # implies that L*v should be exact for monomials up to order 3.
+            @testset "4th order" begin
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+			    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+                Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+                # 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
+                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
+                @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
+                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
+            end
+        end
+
+        @testset "2D" begin
+            l2(v) = sqrt(prod(spacing(g_2D))*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))
+            end
+            v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
+            v_yy = evalOn(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.
+            @testset "2nd order" begin
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+                closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+                Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,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*v ≈ v_yy rtol = 5e-2 norm = l2
+            end
+
+            # 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
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+                closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+                Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,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
+                @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*v ≈ v_yy rtol = 5e-4 norm = l2
+            end
+        end
+    end
+end
--- a/test/SbpOperators/volumeops/derivatives/secondderivative_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,108 +0,0 @@
-using Test
-
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
-
-import Sbplib.SbpOperators.VolumeOperator
-
-@testset "SecondDerivative" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    Lx = 3.5
-    Ly = 3.
-    g_1D = EquidistantGrid(121, 0.0, Lx)
-    g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
-
-    @testset "Constructors" begin
-        @testset "1D" begin
-            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ₓₓ = 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
-        end
-    end
-
-    # Exact differentiation is measured point-wise. In other cases
-    # the error is measured in the l2-norm.
-    @testset "Accuracy" begin
-        @testset "1D" begin
-            l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
-            monomials = ()
-            maxOrder = 4;
-            for i = 0:maxOrder-1
-                f_i(x) = 1/factorial(i)*x^i
-                monomials = (monomials...,evalOn(g_1D,f_i))
-            end
-            v = evalOn(g_1D,x -> sin(x))
-            vₓₓ = evalOn(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.
-            @testset "2nd order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                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
-                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
-            end
-
-            # 4th order interior stencil, 2nd order boundary stencil,
-            # 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ₓₓ = 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
-                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
-                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
-                @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
-                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
-            end
-        end
-
-        @testset "2D" begin
-            l2(v) = sqrt(prod(spacing(g_2D))*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))
-            end
-            v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
-            v_yy = evalOn(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.
-            @testset "2nd order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=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
-                @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
-            end
-
-            # 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
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                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
-                @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*v ≈ v_yy rtol = 5e-4 norm = l2
-            end
-        end
-    end
-end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -14,36 +14,39 @@
     g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
     integral(H,v) = sum(H*v)
     @testset "inner_product" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        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}(),op.quadratureClosure)
+            H = inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
             @test H == IdentityMapping{Float64}()
             @test H isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
-            H = inner_product(g_1D,op.quadratureClosure)
-            inner_stencil = CenteredStencil(1.)
-            @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
+            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+            @test H == inner_product(g_1D, quadrature_interior, quadrature_closure)
             @test H isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            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)
+            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 == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
     end
 
     @testset "Sizes" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        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,op.quadratureClosure)
+            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = inner_product(g_2D,op.quadratureClosure)
+            H = inner_product(g_2D, quadrature_interior, quadrature_closure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -59,8 +62,10 @@
             u = evalOn(g_1D,x->sin(x))
 
             @testset "2nd order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_1D,op.quadratureClosure)
+                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)
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -68,8 +73,10 @@
             end
 
             @testset "4th order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = inner_product(g_1D,op.quadratureClosure)
+                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)
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -82,14 +89,18 @@
             v = b*ones(Float64, size(g_2D))
             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 = inner_product(g_2D,op.quadratureClosure)
+                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
             end
             @testset "4th order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = inner_product(g_2D,op.quadratureClosure)
+                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
             end
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -12,40 +12,38 @@
     g_1D = EquidistantGrid(77, 0.0, Lx)
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
     @testset "inverse_inner_product" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        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}(),op.quadratureClosure)
+            Hi = inverse_inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
             @test Hi == IdentityMapping{Float64}()
             @test Hi isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
-            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 == inverse_inner_product(g_1D,closures,inner_stencil)
+            Hi = inverse_inner_product(g_1D,  quadrature_interior, quadrature_closure)
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            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)
+            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 == Hi_x⊗Hi_y
             @test Hi isa TensorMapping{T,2,2} where T
         end
     end
 
     @testset "Sizes" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        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,op.quadratureClosure)
+            Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
             @test domain_size(Hi) == size(g_1D)
             @test range_size(Hi) == size(g_1D)
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -56,16 +54,20 @@
             v = evalOn(g_1D,x->sin(x))
             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 = inner_product(g_1D,op.quadratureClosure)
-                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+                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)
                 @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 = inner_product(g_1D,op.quadratureClosure)
-                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+                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)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -74,16 +76,20 @@
             v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             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 = inner_product(g_2D,op.quadratureClosure)
-                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+                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)
                 @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 = inner_product(g_2D,op.quadratureClosure)
-                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+                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)
                 @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	Thu Jan 20 21:51:53 2022 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -6,169 +6,188 @@
 using Sbplib.RegionIndices
 using Sbplib.StaticDicts
 
+operator_path = sbp_operators_path()*"standard_diagonal.toml"
+# Default stencils (4th order)
+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"])
+e_closure = parse_stencil(stencil_set["e"]["closure"])
+d_closure = parse_stencil(stencil_set["d1"]["closure"])
+quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
+quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
+
 @testset "Laplace" begin
     g_1D = EquidistantGrid(101, 0.0, 1.)
     g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     @testset "Constructors" begin
 
         @testset "1D" begin
-            # Create all tensor mappings included in Laplace
-            Δ = laplace(g_1D, op.innerStencil, op.closureStencils)
-            H = inner_product(g_1D, op.quadratureClosure)
-            Hi = inverse_inner_product(g_1D, op.quadratureClosure)
+            Δ = laplace(g_1D, inner_stencil, closure_stencils)
+            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+            Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
 
             (id_l, id_r) = boundary_identifiers(g_1D)
 
-            e_l = boundary_restriction(g_1D,op.eClosure,id_l)
-            e_r = boundary_restriction(g_1D,op.eClosure,id_r)
+            e_l = boundary_restriction(g_1D, e_closure,id_l)
+            e_r = boundary_restriction(g_1D, e_closure,id_r)
             e_dict = StaticDict(id_l => e_l, id_r => e_r)
 
-            d_l = normal_derivative(g_1D,op.dClosure,id_l)
-            d_r = normal_derivative(g_1D,op.dClosure,id_r)
+            d_l = normal_derivative(g_1D, d_closure,id_l)
+            d_r = normal_derivative(g_1D, d_closure,id_r)
             d_dict = StaticDict(id_l => d_l, id_r => d_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)
+            H_l = inner_product(boundary_grid(g_1D,id_l), quadrature_interior, quadrature_closure)
+            H_r = inner_product(boundary_grid(g_1D,id_r), quadrature_interior, quadrature_closure)
             Hb_dict = StaticDict(id_l => H_l, id_r => H_r)
 
-            L = Laplace(g_1D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            @test L == Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
+            L = Laplace(g_1D, operator_path; order=4)
+            @test L == Laplace(Δ, H, Hi, e_dict, d_dict, Hb_dict)
             @test L isa TensorMapping{T,1,1}  where T
-            @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
+            @inferred Laplace(Δ, H, Hi, e_dict, d_dict, Hb_dict)
+            # REVIEW: The tests above seem very tied to the implementation. Is
+            # it important that the components of the operator set are stored
+            # in static dicts? Is something like below better?
+            #
+            # ```
+            # L = Laplace(g_1D, operator_path; order=4)
+            # @test L isa TensorMapping{T,1,1}  where T
+            # @test boundary_restriction(L,id_l) ==  boundary_restriction(g_1D, e_closure,id_l)
+            # ...
+            # ```
+            # I guess this is more or less simply a reorganization of the test and skipping testing for the struct layout
         end
         @testset "3D" begin
-            # Create all tensor mappings included in Laplace
-            Δ = laplace(g_3D, op.innerStencil, op.closureStencils)
-            H = inner_product(g_3D, op.quadratureClosure)
-            Hi = inverse_inner_product(g_3D, op.quadratureClosure)
+            Δ = laplace(g_3D, inner_stencil, closure_stencils)
+            H = inner_product(g_3D, quadrature_interior, quadrature_closure)
+            Hi = inverse_inner_product(g_3D, quadrature_interior, quadrature_closure)
 
             (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
-            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_l = boundary_restriction(g_3D, e_closure,id_l)
+            e_r = boundary_restriction(g_3D, e_closure,id_r)
+            e_s = boundary_restriction(g_3D, e_closure,id_s)
+            e_n = boundary_restriction(g_3D, e_closure,id_n)
+            e_b = boundary_restriction(g_3D, e_closure,id_b)
+            e_t = boundary_restriction(g_3D, e_closure,id_t)
             e_dict = StaticDict(id_l => e_l, id_r => e_r,
                                 id_s => e_s, id_n => e_n,
                                 id_b => e_b, id_t => e_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_l = normal_derivative(g_3D, d_closure,id_l)
+            d_r = normal_derivative(g_3D, d_closure,id_r)
+            d_s = normal_derivative(g_3D, d_closure,id_s)
+            d_n = normal_derivative(g_3D, d_closure,id_n)
+            d_b = normal_derivative(g_3D, d_closure,id_b)
+            d_t = normal_derivative(g_3D, d_closure,id_t)
             d_dict = StaticDict(id_l => d_l, id_r => d_r,
                                 id_s => d_s, id_n => d_n,
                                 id_b => d_b, id_t => d_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)
+            H_l = inner_product(boundary_grid(g_3D,id_l), quadrature_interior, quadrature_closure)
+            H_r = inner_product(boundary_grid(g_3D,id_r), quadrature_interior, quadrature_closure)
+            H_s = inner_product(boundary_grid(g_3D,id_s), quadrature_interior, quadrature_closure)
+            H_n = inner_product(boundary_grid(g_3D,id_n), quadrature_interior, quadrature_closure)
+            H_b = inner_product(boundary_grid(g_3D,id_b), quadrature_interior, quadrature_closure)
+            H_t = inner_product(boundary_grid(g_3D,id_t), quadrature_interior, quadrature_closure)
             Hb_dict = StaticDict(id_l => H_l, id_r => H_r,
                                  id_s => H_s, id_n => H_n,
                                  id_b => H_b, id_t => H_t)
 
-            L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            L = Laplace(g_3D, operator_path; order=4)
             @test L == Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
             @test L isa TensorMapping{T,3,3} where T
             @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
         end
     end
 
+    # REVIEW: Is this testset misplaced? Should it really be inside the "Laplace" testset?
     @testset "laplace" begin
         @testset "1D" begin
-            L = laplace(g_1D, op.innerStencil, op.closureStencils)
-            @test L == second_derivative(g_1D, op.innerStencil, op.closureStencils)
+            L = laplace(g_1D, inner_stencil, closure_stencils)
+            @test L == second_derivative(g_1D, inner_stencil, closure_stencils)
             @test L isa TensorMapping{T,1,1}  where T
         end
         @testset "3D" begin
-            L = laplace(g_3D, op.innerStencil, op.closureStencils)
+            L = laplace(g_3D, inner_stencil, closure_stencils)
             @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)
+            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)
             @test L == Dxx + Dyy + Dzz
             @test L isa TensorMapping{T,3,3} where T
         end
     end
 
     @testset "inner_product" begin
-        L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        @test inner_product(L) == inner_product(g_3D,op.quadratureClosure)
+        L = Laplace(g_3D, operator_path; order=4)
+        @test inner_product(L) == inner_product(g_3D, quadrature_interior, quadrature_closure)
     end
 
     @testset "inverse_inner_product" begin
-        L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        @test inverse_inner_product(L) == inverse_inner_product(g_3D,op.quadratureClosure)
+        L = Laplace(g_3D, operator_path; order=4)
+        @test inverse_inner_product(L) == inverse_inner_product(g_3D, quadrature_interior, quadrature_closure)
     end
 
     @testset "boundary_restriction" begin
-        L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        L = Laplace(g_3D, operator_path; order=4)
         (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
-        @test boundary_restriction(L,id_l) == boundary_restriction(g_3D,op.eClosure,id_l)
-        @test boundary_restriction(L,id_r) == boundary_restriction(g_3D,op.eClosure,id_r)
-        @test boundary_restriction(L,id_s) == boundary_restriction(g_3D,op.eClosure,id_s)
-        @test boundary_restriction(L,id_n) == boundary_restriction(g_3D,op.eClosure,id_n)
-        @test boundary_restriction(L,id_b) == boundary_restriction(g_3D,op.eClosure,id_b)
-        @test boundary_restriction(L,id_t) == boundary_restriction(g_3D,op.eClosure,id_t)
+        @test boundary_restriction(L, id_l) == boundary_restriction(g_3D, e_closure,id_l)
+        @test boundary_restriction(L, id_r) == boundary_restriction(g_3D, e_closure,id_r)
+        @test boundary_restriction(L, id_s) == boundary_restriction(g_3D, e_closure,id_s)
+        @test boundary_restriction(L, id_n) == boundary_restriction(g_3D, e_closure,id_n)
+        @test boundary_restriction(L, id_b) == boundary_restriction(g_3D, e_closure,id_b)
+        @test boundary_restriction(L, id_t) == boundary_restriction(g_3D, e_closure,id_t)
 
         ids = boundary_identifiers(g_3D)
-        es = boundary_restriction(L,ids)
-        @test es ==  (boundary_restriction(L,id_l),
-                      boundary_restriction(L,id_r),
-                      boundary_restriction(L,id_s),
-                      boundary_restriction(L,id_n),
-                      boundary_restriction(L,id_b),
-                      boundary_restriction(L,id_t));
-        @test es == boundary_restriction(L,ids...)
+        es = boundary_restriction(L, ids)
+        @test es ==  (boundary_restriction(L, id_l),
+                      boundary_restriction(L, id_r),
+                      boundary_restriction(L, id_s),
+                      boundary_restriction(L, id_n),
+                      boundary_restriction(L, id_b),
+                      boundary_restriction(L, id_t));
+        @test es == boundary_restriction(L, ids...)
     end
 
     @testset "normal_derivative" begin
-        L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        L = Laplace(g_3D, operator_path; order=4)
         (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
-        @test normal_derivative(L,id_l) == normal_derivative(g_3D,op.dClosure,id_l)
-        @test normal_derivative(L,id_r) == normal_derivative(g_3D,op.dClosure,id_r)
-        @test normal_derivative(L,id_s) == normal_derivative(g_3D,op.dClosure,id_s)
-        @test normal_derivative(L,id_n) == normal_derivative(g_3D,op.dClosure,id_n)
-        @test normal_derivative(L,id_b) == normal_derivative(g_3D,op.dClosure,id_b)
-        @test normal_derivative(L,id_t) == normal_derivative(g_3D,op.dClosure,id_t)
+        @test normal_derivative(L, id_l) == normal_derivative(g_3D, d_closure,id_l)
+        @test normal_derivative(L, id_r) == normal_derivative(g_3D, d_closure,id_r)
+        @test normal_derivative(L, id_s) == normal_derivative(g_3D, d_closure,id_s)
+        @test normal_derivative(L, id_n) == normal_derivative(g_3D, d_closure,id_n)
+        @test normal_derivative(L, id_b) == normal_derivative(g_3D, d_closure,id_b)
+        @test normal_derivative(L, id_t) == normal_derivative(g_3D, d_closure,id_t)
 
         ids = boundary_identifiers(g_3D)
-        ds = normal_derivative(L,ids)
-        @test ds ==  (normal_derivative(L,id_l),
-                      normal_derivative(L,id_r),
-                      normal_derivative(L,id_s),
-                      normal_derivative(L,id_n),
-                      normal_derivative(L,id_b),
-                      normal_derivative(L,id_t));
-        @test ds == normal_derivative(L,ids...)
+        ds = normal_derivative(L, ids)
+        @test ds ==  (normal_derivative(L, id_l),
+                      normal_derivative(L, id_r),
+                      normal_derivative(L, id_s),
+                      normal_derivative(L, id_n),
+                      normal_derivative(L, id_b),
+                      normal_derivative(L, id_t));
+        @test ds == normal_derivative(L, ids...)
     end
 
     @testset "boundary_quadrature" begin
-        L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        L = Laplace(g_3D, operator_path; order=4)
         (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
-        @test boundary_quadrature(L,id_l) == inner_product(boundary_grid(g_3D,id_l),op.quadratureClosure)
-        @test boundary_quadrature(L,id_r) == inner_product(boundary_grid(g_3D,id_r),op.quadratureClosure)
-        @test boundary_quadrature(L,id_s) == inner_product(boundary_grid(g_3D,id_s),op.quadratureClosure)
-        @test boundary_quadrature(L,id_n) == inner_product(boundary_grid(g_3D,id_n),op.quadratureClosure)
-        @test boundary_quadrature(L,id_b) == inner_product(boundary_grid(g_3D,id_b),op.quadratureClosure)
-        @test boundary_quadrature(L,id_t) == inner_product(boundary_grid(g_3D,id_t),op.quadratureClosure)
+        @test boundary_quadrature(L, id_l) == inner_product(boundary_grid(g_3D, id_l), quadrature_interior, quadrature_closure)
+        @test boundary_quadrature(L, id_r) == inner_product(boundary_grid(g_3D, id_r), quadrature_interior, quadrature_closure)
+        @test boundary_quadrature(L, id_s) == inner_product(boundary_grid(g_3D, id_s), quadrature_interior, quadrature_closure)
+        @test boundary_quadrature(L, id_n) == inner_product(boundary_grid(g_3D, id_n), quadrature_interior, quadrature_closure)
+        @test boundary_quadrature(L, id_b) == inner_product(boundary_grid(g_3D, id_b), quadrature_interior, quadrature_closure)
+        @test boundary_quadrature(L, id_t) == inner_product(boundary_grid(g_3D, id_t), quadrature_interior, quadrature_closure)
 
         ids = boundary_identifiers(g_3D)
-        H_gammas = boundary_quadrature(L,ids)
-        @test H_gammas ==  (boundary_quadrature(L,id_l),
-                            boundary_quadrature(L,id_r),
-                            boundary_quadrature(L,id_s),
-                            boundary_quadrature(L,id_n),
-                            boundary_quadrature(L,id_b),
-                            boundary_quadrature(L,id_t));
-        @test H_gammas == boundary_quadrature(L,ids...)
+        H_gammas = boundary_quadrature(L, ids)
+        @test H_gammas ==  (boundary_quadrature(L, id_l),
+                            boundary_quadrature(L, id_r),
+                            boundary_quadrature(L, id_s),
+                            boundary_quadrature(L, id_n),
+                            boundary_quadrature(L, id_b),
+                            boundary_quadrature(L, id_t));
+        @test H_gammas == boundary_quadrature(L, ids...)
     end
 
     # Exact differentiation is measured point-wise. In other cases
@@ -187,7 +206,10 @@
         # 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)
+            stencil_set = read_stencil_set(operator_path; order=2)
+            inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+            closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+            L = laplace(g_3D, inner_stencil, closure_stencils)
             @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
@@ -197,7 +219,10 @@
         # 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)
+            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"])
+            L = laplace(g_3D, inner_stencil, closure_stencils)
             # 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