diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index f925133..e3ee894 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1,6 +1,5 @@
 # See https://domluna.github.io/JuliaFormatter.jl/stable/ for a list of options
 style = "blue"
-
 margin = 92
 
 align_assignment = true
@@ -15,4 +14,4 @@ format_docstrings = true
 conditional_to_if = true
 indent_submodule = true
 
-always_for_in=true
\ No newline at end of file
+always_for_in = true
diff --git a/.github/codecov.yml b/.github/codecov.yml
new file mode 100644
index 0000000..d68eef7
--- /dev/null
+++ b/.github/codecov.yml
@@ -0,0 +1,12 @@
+coverage:
+  status:
+    project:
+      default:
+        informational: true
+    patch:
+      default:
+        informational: true
+
+comment:
+  layout: "reach, diff, flags, files"
+  require_changes: true
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 0000000..f207934
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,12 @@
+version: 2
+
+updates:
+  - package-ecosystem: "julia"
+    directory: "/"
+    schedule:
+      interval: "weekly"
+
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      interval: "weekly"
diff --git a/.github/workflows/AutoMerge.yml b/.github/workflows/AutoMerge.yml
new file mode 100644
index 0000000..2fb156a
--- /dev/null
+++ b/.github/workflows/AutoMerge.yml
@@ -0,0 +1,23 @@
+name: Auto-merge dependency PRs
+
+on: pull_request_target
+
+permissions:
+  contents: write
+  pull-requests: write
+
+jobs:
+  dependabot:
+    name: Dependabot
+    runs-on: ubuntu-latest
+    if: github.actor == 'dependabot[bot]'
+    steps:
+      - name: Fetch Dependabot metadata
+        uses: dependabot/fetch-metadata@v3
+        with:
+          github-token: ${{ secrets.GITHUB_TOKEN }}
+      - name: Enable auto-merge
+        run: gh pr merge --auto --squash "$PR_URL"
+        env:
+          PR_URL: ${{ github.event.pull_request.html_url }}
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 522f9a0..2ffa0d1 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -1,41 +1,64 @@
 name: CI
+
+env:
+  JULIA_PKG_SERVER_REGISTRY_PREFERENCE: conservative
+
 on:
   push:
     branches:
-      - 'master'
-      - 'dev'
-    tags: '*'
+      - main
+      - dev
+    tags: "*"
   pull_request:
     branches:
-      - 'master'
+      - main
+      - dev
   release:
     types: [published]
 
+permissions:
+  actions: write
+  contents: read
+
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
         version:
-          - '1.10'
+          - "1.10"
+          - "1"
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v3
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v3
+      - name: Ensure General registry
+        run: |
+          julia -e '
+            using Pkg
+            if isempty(Pkg.Registry.reachable_registries())
+              try
+                Pkg.Registry.add()
+              catch
+                ENV["JULIA_PKG_SERVER"] = ""
+                Pkg.Registry.add(Pkg.RegistrySpec(
+                  url="https://github.com/JuliaRegistries/General.git",
+                ))
+              end
+            end'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v6
         with:
           files: lcov.info
-        env:
-          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
\ No newline at end of file
+          token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index d66415e..79c907f 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -1,30 +1,49 @@
 name: Documentation
 
+env:
+  JULIA_PKG_SERVER_REGISTRY_PREFERENCE: conservative
+
 on:
   push:
     branches:
-      - 'master'
-      - 'dev'
-    tags: '*'
+      - main
+      - dev
+    tags: "*"
   pull_request:
   release:
     types: [published]
 
 jobs:
   build:
+    name: Documentation
     permissions:
       contents: write
       statuses: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v3
         with:
-          version: '1'
+          version: "1"
+      - uses: julia-actions/cache@v3
+      - name: Ensure General registry
+        run: |
+          julia -e '
+            using Pkg
+            if isempty(Pkg.Registry.reachable_registries())
+              try
+                Pkg.Registry.add()
+              catch
+                ENV["JULIA_PKG_SERVER"] = ""
+                Pkg.Registry.add(Pkg.RegistrySpec(
+                  url="https://github.com/JuliaRegistries/General.git",
+                ))
+              end
+            end'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
         env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
\ No newline at end of file
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+        run: julia --project=docs/ docs/make.jl
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
index f49313b..ef6819f 100644
--- a/.github/workflows/TagBot.yml
+++ b/.github/workflows/TagBot.yml
@@ -1,9 +1,14 @@
 name: TagBot
+
 on:
   issue_comment:
     types:
       - created
   workflow_dispatch:
+
+permissions:
+  contents: write
+
 jobs:
   TagBot:
     if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
diff --git a/.gitignore b/.gitignore
index 3b5673d..6b6bc50 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,8 +1,12 @@
+.vscode/
+*.dat
+*.jld2
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
+*.mem
+lcov.info
 /Manifest.toml
 /docs/Manifest.toml
 /docs/build/
-lcov.info
-*.vscode
\ No newline at end of file
+/docs/site/
diff --git a/CITATION.cff b/CITATION.cff
new file mode 100644
index 0000000..008d813
--- /dev/null
+++ b/CITATION.cff
@@ -0,0 +1,12 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it using the metadata below."
+title: "CorrectionFactorMatrixMethod.jl"
+type: software
+authors:
+  - family-names: Tetzner
+    given-names: Joshua
+repository-code: "https://github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl"
+url: "https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl/"
+license: MIT
+version: 0.1.0
+# doi: "Add the Zenodo DOI after the first archived release."
diff --git a/Project.toml b/Project.toml
index d23f2af..3f131dc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "CorrectionFactorMatrixMethod"
 uuid = "1b32fdd5-07c8-465a-bf7d-415ac82a78be"
 version = "0.1.0"
-authors = ["""Joshua <joshua.tetzner@uni-rostock.de\n>"""]
+authors = ["Joshua Tetzner <joshua.tetzner@uni-rostock.de>"]
 
 [deps]
 BlockSparseMatrices = "1ef78180-cc57-4e70-b6ec-c3e471307147"
@@ -22,8 +22,47 @@ CFMMBEAST = ["BEAST", "CompScienceMeshes"]
 CFMMExaFMMt = ["ExaFMMt", "BEAST"]
 
 [compat]
-BlockSparseMatrices = "0.3.0"
-CompScienceMeshes = "0.12"
-MKL = "0.6.3"
-OhMyThreads = "0.8.4"
-SparseArrays = "1.12.0"
+Aqua = "0.8"
+BEAST = "2"
+BlockSparseMatrices = "0.3"
+CompScienceMeshes = "0.11, 0.12"
+ExaFMMt = "0.1"
+ExplicitImports = "1"
+H2Trees = "0.3"
+JET = "0.9, 0.10, 0.11"
+JuliaFormatter = "2"
+LinearAlgebra = "1.10"
+LinearMaps = "3"
+MKL = "0.9"
+OhMyThreads = "0.8"
+SparseArrays = "1.10"
+StaticArrays = "1"
+Test = "1.10"
+TestItemRunner = "1"
+TestItems = "1"
+julia = "1.10"
+
+[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
+
+[targets]
+test = [
+    "Aqua",
+    "BEAST",
+    "CompScienceMeshes",
+    "ExaFMMt",
+    "ExplicitImports",
+    "JET",
+    "JuliaFormatter",
+    "StaticArrays",
+    "Test",
+    "TestItemRunner",
+    "TestItems",
+]
diff --git a/README.md b/README.md
index 3786126..37bad6e 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,80 @@
 # CorrectionFactorMatrixMethod.jl
+
+[![Docs-stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl/stable/)
+[![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl/dev/)
+[![MIT license](https://img.shields.io/badge/License-MIT-blue.svg)](https://github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl/blob/main/LICENSE)
+[![CI](https://github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl/actions/workflows/CI.yml)
+[![codecov](https://codecov.io/gh/JoshuaTetzner/CorrectionFactorMatrixMethod.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JoshuaTetzner/CorrectionFactorMatrixMethod.jl)
+
+## Introduction
+
+This package builds matrix-free boundary-element operators that combine a fast
+multipole approximation of the far interactions with directly assembled
+near-field corrections. The result behaves like the dense boundary-element
+matrix under matrix-vector products, but is never assembled explicitly, so it
+can be used directly inside an iterative solver such as `Krylov.gmres`.
+
+The generic correction and linear-map machinery lives in the core package.
+Support for concrete operators and spaces is provided through package
+extensions for
+[BEAST.jl](https://github.com/krcools/BEAST.jl),
+[CompScienceMeshes.jl](https://github.com/krcools/CompScienceMeshes.jl), and
+[ExaFMMt.jl](https://github.com/JoshuaTetzner/ExaFMMt.jl), which are loaded
+automatically once these packages are available.
+
+## Correction-Factor Matrix Method
+
+For a boundary-element matrix $A$, the far field is approximated by an
+FMM-backed map $A_\mathrm{FMM}$. The near interactions are evaluated with the
+boundary-element quadrature and corrected for the part already represented by
+the FMM:
+
+$$Ax \approx A_\mathrm{FMM}\,x + \left(A_\mathrm{near} - A_\mathrm{FMM,near}\right)x.$$
+
+The sparse correction blocks are selected from an
+[H2Trees.jl](https://github.com/djukic14/H2Trees.jl) tree. A matrix-vector
+product evaluates the FMM map first and then adds the sparse near correction.
+Transpose and adjoint products follow the corresponding `LinearMaps.jl`
+interfaces. Further details are given in the
+[documentation](https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl/dev/details/method/).
+
+## Installation
+
+Installing CorrectionFactorMatrixMethod is done by entering the package manager
+(enter `]` at the Julia REPL) and issuing:
+
+```
+pkg> add https://github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl.git
+```
+
+The supported boundary-element integration requires Julia 1.10 or later.
+
+## First steps
+
+Load the optional BEAST, CompScienceMeshes, and ExaFMMt packages, build a
+boundary-element operator with its test and trial spaces, and assemble the
+correction-factor operator:
+
+```julia
+using BEAST, CompScienceMeshes, ExaFMMt
+using CorrectionFactorMatrixMethod
+
+mesh = meshsphere(1.0, 0.4)
+space = raviartthomas(mesh)
+operator = Maxwell3D.singlelayer(; wavenumber=1.0)
+
+matrix = CFMM.assemble(operator, space)
+result = matrix * rand(scalartype(operator), numfunctions(space))
+```
+
+[`CFMM.assemble`](https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl/dev/manual/manual/)
+constructs an optimized tree automatically. Existing trees and all lower-level
+FMM, quadrature, and scheduler options can be supplied as keywords. The
+returned operator implements the `LinearMaps.jl` interface, so it can be passed
+straight to an iterative solver. Runnable EFIE and MFIE examples are in the
+[`examples/`](examples) directory and in the
+[documentation](https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl/dev/).
+
+## References
+
+- [1] Adelman, Ross, Nail A. Gumerov, and Ramani Duraiswami. *FMM/GPU-Accelerated Boundary Element Method for Computational Magnetics and Electrostatics.* IEEE Transactions on Magnetics 53, no. 12 (December 2017): 1–11. [https://doi.org/10.1109/TMAG.2017.2725951](https://doi.org/10.1109/TMAG.2017.2725951).
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..59e9900
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+CorrectionFactorMatrixMethod = "1b32fdd5-07c8-465a-bf7d-415ac82a78be"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+
+[compat]
+CorrectionFactorMatrixMethod = "0.1"
+Documenter = "1"
+julia = "1.10"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..cd0608f
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,45 @@
+using CorrectionFactorMatrixMethod
+using Documenter
+
+DocMeta.setdocmeta!(
+    CorrectionFactorMatrixMethod,
+    :DocTestSetup,
+    :(using CorrectionFactorMatrixMethod);
+    recursive=true,
+)
+
+makedocs(;
+    modules=[CorrectionFactorMatrixMethod, CorrectionFactorMatrixMethod.CFMM],
+    authors="Joshua Tetzner",
+    repo=Remotes.GitHub("JoshuaTetzner", "CorrectionFactorMatrixMethod.jl"),
+    sitename="CorrectionFactorMatrixMethod.jl",
+    format=Documenter.HTML(;
+        canonical="https://JoshuaTetzner.github.io/CorrectionFactorMatrixMethod.jl",
+        edit_link="dev",
+        assets=String[],
+    ),
+    pages=[
+        "Introduction" => "index.md",
+        "Manual" =>
+            ["General Usage" => "manual/manual.md", "Examples" => "manual/examples.md"],
+        "Details" => [
+            "Correction-Factor Method" => "details/method.md",
+            "Package Extensions" => "details/extensions.md",
+        ],
+        "Contributing" => "contributing.md",
+        "API Reference" => "apiref.md",
+    ],
+    checkdocs=:exports,
+)
+
+for (root, _, files) in walkdir(joinpath(@__DIR__, "build"))
+    for file in files
+        endswith(file, ".dat") && rm(joinpath(root, file))
+    end
+end
+
+deploydocs(;
+    repo="github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl.git",
+    devbranch="dev",
+    push_preview=false,
+)
diff --git a/docs/src/apiref.md b/docs/src/apiref.md
new file mode 100644
index 0000000..4740c28
--- /dev/null
+++ b/docs/src/apiref.md
@@ -0,0 +1,12 @@
+# API Reference
+
+```@docs
+CorrectionFactorMatrixMethod
+CorrectionFactorMatrixMethod.CFMM
+```
+
+```@autodocs
+Modules = [CorrectionFactorMatrixMethod, CorrectionFactorMatrixMethod.CFMM]
+Private = true
+Order = [:type, :function]
+```
diff --git a/docs/src/contributing.md b/docs/src/contributing.md
new file mode 100644
index 0000000..9f0c262
--- /dev/null
+++ b/docs/src/contributing.md
@@ -0,0 +1,23 @@
+# Contributing
+
+Development work targets the `dev` branch. Open a pull request from `dev` to
+`main` for stable changes, and merge it with a merge commit.
+
+Run the package checks locally with:
+
+```julia
+using Pkg
+Pkg.test()
+```
+
+Build the documentation with:
+
+```julia
+using Pkg
+Pkg.develop(PackageSpec(path=pwd()))
+Pkg.instantiate()
+include("docs/make.jl")
+```
+
+Code must pass the formatter, Aqua, JET, and ExplicitImports checks included in
+the test suite.
diff --git a/docs/src/details/extensions.md b/docs/src/details/extensions.md
new file mode 100644
index 0000000..6504471
--- /dev/null
+++ b/docs/src/details/extensions.md
@@ -0,0 +1,20 @@
+# Package Extensions
+
+CorrectionFactorMatrixMethod.jl uses package extensions to keep optional
+boundary-element and FMM dependencies out of the core load path.
+
+`CFMMBEAST` loads when `BEAST` and `CompScienceMeshes` are available. It
+provides quadrature, collocation, corrected-kernel assembly, and supported
+Helmholtz and Maxwell operator methods.
+
+`CFMMExaFMMt` loads when `ExaFMMt` and `BEAST` are available. It maps BEAST
+operator parameters to ExaFMMt options and constructs the FMM maps used by
+[`CFMM.assemble`](@ref).
+
+Load all three optional packages before constructing a supported
+boundary-element operator:
+
+```julia
+using BEAST, CompScienceMeshes, ExaFMMt
+using CorrectionFactorMatrixMethod
+```
diff --git a/docs/src/details/method.md b/docs/src/details/method.md
new file mode 100644
index 0000000..38dfc24
--- /dev/null
+++ b/docs/src/details/method.md
@@ -0,0 +1,16 @@
+# Correction-Factor Method
+
+For a boundary-element matrix ``A``, the package approximates the far field by
+an FMM-backed map ``A_{\mathrm{FMM}}``. Near interactions are evaluated with
+the boundary-element quadrature and corrected for the interactions already
+represented by the FMM:
+
+```math
+A x \approx A_{\mathrm{FMM}}x +
+\left(A_{\mathrm{near}} - A_{\mathrm{FMM,near}}\right)x.
+```
+
+The sparse correction blocks are selected from an `H2Trees.BlockTree`.
+Matrix-vector products evaluate the FMM map first and then add the sparse near
+correction. Transpose and adjoint products follow the corresponding
+`LinearMaps.jl` interfaces.
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..438fe4b
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,21 @@
+# CorrectionFactorMatrixMethod.jl
+
+CorrectionFactorMatrixMethod.jl builds matrix-free boundary-element operators
+that combine a fast multipole approximation of far interactions with directly
+assembled near-field corrections.
+
+The package keeps the generic correction and linear-map machinery in its core.
+Optional integrations with
+[BEAST.jl](https://github.com/krcools/BEAST.jl),
+[CompScienceMeshes.jl](https://github.com/krcools/CompScienceMeshes.jl), and
+[ExaFMMt.jl](https://github.com/JoshuaTetzner/ExaFMMt.jl) are loaded through
+Julia package extensions.
+
+## Installation
+
+```julia
+using Pkg
+Pkg.add(url="https://github.com/JoshuaTetzner/CorrectionFactorMatrixMethod.jl")
+```
+
+The current integration requires Julia 1.10 or later.
diff --git a/docs/src/manual/examples.md b/docs/src/manual/examples.md
new file mode 100644
index 0000000..9b84c2b
--- /dev/null
+++ b/docs/src/manual/examples.md
@@ -0,0 +1,84 @@
+# Examples
+
+The examples use BEAST for the boundary-element spaces and operators, ExaFMMt
+for the far interactions, and `Krylov.gmres` for the matrix-free solve. Install
+the example environment and run either script from the package root:
+
+```sh
+julia --project=examples -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+julia --project=examples examples/efie.jl
+julia --project=examples examples/mfie.jl
+```
+
+Both scripts write an interactive HTML plot containing the far-field pattern,
+the electric field in the ``yz`` plane, and the surface-current magnitude.
+They use the following packages:
+
+```julia
+using BEAST
+using CompScienceMeshes
+using CorrectionFactorMatrixMethod
+using ExaFMMt
+using Krylov
+using LinearAlgebra
+using PlotlyJS
+```
+
+## EFIE
+
+The electric-field integral equation uses the Maxwell single-layer operator
+with the same Raviart--Thomas space for testing and trial functions:
+
+```julia
+mesh = meshsphere(1.0, 0.4)
+space = raviartthomas(mesh)
+
+wavenumber = 1.0
+operator = Maxwell3D.singlelayer(; wavenumber)
+excitation = Maxwell3D.planewave(;
+    direction=ẑ, polarization=x̂, wavenumber
+)
+rhs = assemble((n × excitation) × n, space)
+
+matrix = CFMM.assemble(operator, space)
+current, stats = Krylov.gmres(
+    matrix, rhs; rtol=1.0e-4, itmax=200, history=true, verbose=1
+)
+```
+
+The complete runnable example, including the plots, is in
+`examples/efie.jl`.
+
+## MFIE
+
+For the magnetic-field integral equation, the nonlocal double-layer operator
+uses CFMM while the local ``\frac{1}{2}I`` contribution is assembled directly.
+Buffa--Christiansen functions provide the test space:
+
+```julia
+mesh = meshsphere(1.0, 0.4)
+trialspace = raviartthomas(mesh)
+testspace = buffachristiansen(mesh)
+
+permittivity = 1.0
+permeability = 1.0
+frequency = 1.0
+wavenumber = frequency * √(permittivity * permeability)
+
+operator = Maxwell3D.doublelayer(; wavenumber)
+excitation = Maxwell3D.planewave(;
+    direction=ẑ, polarization=x̂, wavenumber
+)
+magneticfield = -1 / (im * permeability * frequency) * curl(excitation)
+rhs = assemble((n × magneticfield) × n, testspace)
+
+doublelayer = CFMM.assemble(operator, testspace, trialspace)
+identity = assemble(NCross(), testspace, trialspace)
+matrix = doublelayer + 0.5 * identity
+current, stats = Krylov.gmres(
+    matrix, rhs; rtol=1.0e-4, itmax=200, history=true, verbose=1
+)
+```
+
+The complete runnable example, including the plots, is in
+`examples/mfie.jl`.
diff --git a/docs/src/manual/manual.md b/docs/src/manual/manual.md
new file mode 100644
index 0000000..b6a5cb3
--- /dev/null
+++ b/docs/src/manual/manual.md
@@ -0,0 +1,26 @@
+# General Usage
+
+The recommended entry point is [`CFMM.assemble`](@ref). It accepts a boundary
+integral operator and its test and trial spaces, constructs an optimized
+`H2Trees.BlockTree`, and returns a [`PetrovGalerkinCFMM`](@ref).
+
+```julia
+matrix = CFMM.assemble(operator, testspace, trialspace)
+```
+
+The lower-level `PetrovGalerkinCFMM` constructor accepts an existing tree.
+Keywords control near and far quadrature, scheduling, and whether a separate
+transpose FMM is built.
+
+The default FMM configuration uses ExaFMMt with expansion order `p = 8` and a
+critical leaf size of `ncrit = 50`.
+
+```@example configuration
+using CorrectionFactorMatrixMethod
+
+configuration = CorrectionFactorMatrixMethod.FMMFunctor(; p=10, ncrit=80)
+(p=configuration.p, ncrit=configuration.ncrit)
+```
+
+The resulting operator implements the `LinearMaps.jl` interface, so it can be
+multiplied by vectors without assembling the dense boundary-element matrix.
diff --git a/examples/Project.toml b/examples/Project.toml
new file mode 100644
index 0000000..91ddfcc
--- /dev/null
+++ b/examples/Project.toml
@@ -0,0 +1,18 @@
+[deps]
+BEAST = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
+CompScienceMeshes = "3e66a162-7b8c-5da0-b8f8-124ecd2c3ae1"
+CorrectionFactorMatrixMethod = "1b32fdd5-07c8-465a-bf7d-415ac82a78be"
+ExaFMMt = "add1291b-d076-47cc-8667-64576de04ee0"
+Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
+
+[compat]
+BEAST = "2"
+CompScienceMeshes = "0.11, 0.12"
+CorrectionFactorMatrixMethod = "0.1"
+ExaFMMt = "0.1"
+Krylov = "0.10"
+LinearAlgebra = "1.10"
+PlotlyJS = "0.18"
+julia = "1.10"
diff --git a/examples/efie.jl b/examples/efie.jl
new file mode 100644
index 0000000..033f476
--- /dev/null
+++ b/examples/efie.jl
@@ -0,0 +1,25 @@
+using BEAST
+using CompScienceMeshes
+using CorrectionFactorMatrixMethod
+using ExaFMMt
+using Krylov
+using LinearAlgebra
+using PlotlyJS
+
+include("plotresults.jl")
+
+mesh = meshsphere(1.0, 0.4)
+space = raviartthomas(mesh)
+
+wavenumber = 1.0
+operator = Maxwell3D.singlelayer(; wavenumber)
+excitation = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber)
+rhs = assemble((n × excitation) × n, space)
+
+matrix = CFMM.assemble(operator, space)
+current, stats = Krylov.gmres(matrix, rhs; rtol=1.0e-4, itmax=200, history=true, verbose=1)
+
+@show stats
+@show norm(matrix * current - rhs) / norm(rhs)
+
+plotresults(current, space, excitation, wavenumber, joinpath(@__DIR__, "efie_results.html"))
diff --git a/examples/mfie.jl b/examples/mfie.jl
new file mode 100644
index 0000000..0a85958
--- /dev/null
+++ b/examples/mfie.jl
@@ -0,0 +1,35 @@
+using BEAST
+using CompScienceMeshes
+using CorrectionFactorMatrixMethod
+using ExaFMMt
+using Krylov
+using LinearAlgebra
+using PlotlyJS
+
+include("plotresults.jl")
+
+mesh = meshsphere(1.0, 0.4)
+trialspace = raviartthomas(mesh)
+testspace = buffachristiansen(mesh)
+
+permittivity = 1.0
+permeability = 1.0
+frequency = 1.0
+wavenumber = frequency * √(permittivity * permeability)
+
+operator = Maxwell3D.doublelayer(; wavenumber)
+excitation = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber)
+magneticfield = -1 / (im * permeability * frequency) * curl(excitation)
+rhs = assemble((n × magneticfield) × n, testspace)
+
+doublelayer = CFMM.assemble(operator, testspace, trialspace)
+identity = assemble(NCross(), testspace, trialspace)
+matrix = doublelayer + 0.5 * identity
+current, stats = Krylov.gmres(matrix, rhs; rtol=1.0e-4, itmax=200, history=true, verbose=1)
+
+@show stats
+@show norm(matrix * current - rhs) / norm(rhs)
+
+plotresults(
+    current, trialspace, excitation, wavenumber, joinpath(@__DIR__, "mfie_results.html")
+)
diff --git a/examples/plotresults.jl b/examples/plotresults.jl
new file mode 100644
index 0000000..2189bb4
--- /dev/null
+++ b/examples/plotresults.jl
@@ -0,0 +1,41 @@
+function plotresults(u, space, excitation, wavenumber, filename)
+    Φ = [0.0]
+    Θ = range(0; stop=π, length=100)
+    points = [point(cos(ϕ) * sin(θ), sin(ϕ) * sin(θ), cos(θ)) for ϕ in Φ for θ in Θ]
+    farfield = potential(MWFarField3D(; wavenumber), points, u, space)
+    currents, geometry = facecurrents(u, space)
+
+    ys = range(-2; stop=2, length=50)
+    zs = range(-4; stop=4, length=100)
+    gridpoints = [point(0, y, z) for y in ys, z in zs]
+    scatteredfield = potential(MWSingleLayerField3D(; wavenumber), gridpoints, u, space)
+    incidentfield = excitation.(gridpoints)
+
+    plot = Plot(
+        Layout(
+            Subplots(;
+                rows=2,
+                cols=2,
+                specs=[Spec() Spec(; rowspan=2); Spec(; kind="mesh3d") missing],
+            ),
+        ),
+    )
+    add_trace!(plot, scatter(; x=Θ, y=norm.(farfield)); row=1, col=1)
+    add_trace!(
+        plot,
+        contour(;
+            x=zs,
+            y=ys,
+            z=norm.(scatteredfield - incidentfield)',
+            colorscale="Viridis",
+            zmin=0,
+            zmax=2,
+            showscale=false,
+        );
+        row=1,
+        col=2,
+    )
+    add_trace!(plot, patch(geometry, norm.(currents); caxis=(0, 2)); row=2, col=1)
+    savefig(plot, filename)
+    return plot
+end
diff --git a/ext/CFMMBEAST/CFMMBEAST.jl b/ext/CFMMBEAST/CFMMBEAST.jl
index 575fec8..2c97b53 100644
--- a/ext/CFMMBEAST/CFMMBEAST.jl
+++ b/ext/CFMMBEAST/CFMMBEAST.jl
@@ -1,10 +1,11 @@
 module CFMMBEAST
 
 using CorrectionFactorMatrixMethod
-import CorrectionFactorMatrixMethod: FMMFunctor, potentials, sources
+import CorrectionFactorMatrixMethod: curlpotentials, divpotentials, normals, potentials
 
 using BEAST
 using CompScienceMeshes
+using H2Trees
 using SparseArrays
 
 include("kernelmatrix.jl")
@@ -35,4 +36,12 @@ function CorrectionFactorMatrixMethod.beta(operator::BEAST.IntegralOperator)
     return operator.beta
 end
 
+function CorrectionFactorMatrixMethod.defaultminhalfsize(
+    testspace::BEAST.Space, trialspace::BEAST.Space
+)
+    _, testhalfsize = H2Trees.boundingbox(BEAST.positions(testspace))
+    _, trialhalfsize = H2Trees.boundingbox(BEAST.positions(trialspace))
+    return max(testhalfsize, trialhalfsize) / 2^10
+end
+
 end
diff --git a/ext/CFMMBEAST/kernelmatrix.jl b/ext/CFMMBEAST/kernelmatrix.jl
index d00cf84..2a6ccec 100644
--- a/ext/CFMMBEAST/kernelmatrix.jl
+++ b/ext/CFMMBEAST/kernelmatrix.jl
@@ -4,7 +4,7 @@ function CorrectionFactorMatrixMethod.AbstractCorrectedKernelMatrix(
     testspace::BEAST.Space,
     trialspace::BEAST.Space;
     nearquadstrat=BEAST.defaultquadstrat(operator, testspace, trialspace),
-    farquadstrat=SafeDoubleNumQStrat(operator, testspace, trialspace),
+    farquadstrat=SafeDoubleNumQStrat(3, 3),
 )
     return CorrectionFactorMatrixMethod.BEASTCorrectedKernelMatrix{scalartype(operator)}(
         BEAST.blockassembler(operator, testspace, trialspace; quadstrat=nearquadstrat),
@@ -31,7 +31,7 @@ function (f::BlockCorrectionFunctor)(v, m, n)
 end
 
 function (blk::CorrectionFactorMatrixMethod.BEASTCorrectedKernelMatrix)(
-    tdata, sdata, matrixblock
+    matrixblock, tdata, sdata
 )
     blk.nearassembler(tdata, sdata, BlockStoreFunctor(matrixblock))
     blk.correctionassembler(tdata, sdata, BlockCorrectionFunctor(matrixblock))
diff --git a/ext/CFMMBEAST/operator/HH3Ddoublelayer.jl b/ext/CFMMBEAST/operator/HH3Ddoublelayer.jl
index 141ed4f..daa6701 100644
--- a/ext/CFMMBEAST/operator/HH3Ddoublelayer.jl
+++ b/ext/CFMMBEAST/operator/HH3Ddoublelayer.jl
@@ -5,11 +5,11 @@ function (fmmfunctor::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     testqp::Matrix,
     trialqp::Matrix,
     fmm;
-    ntasks=Threads.nthreads(),
+    scheduler,
 ) where {T<:BEAST.HH3DDoubleLayerFDBIO}
     Btest, Btrial, normals = HH3DDLpotentialmatrix(testspace, trialspace, testqp, trialqp)
 
-    return CFMMMatrixHH3DDoubleLayer{scalartype(operator)}(
+    return CorrectionFactorMatrixMethod.CFMMMatrixHH3DDoubleLayer{scalartype(operator)}(
         operator, fmm, Btest, Btrial, normals
     )
 end
@@ -17,13 +17,13 @@ end
 function HH3DDLpotentialmatrix(
     testspace::BEAST.Space, trialspace::BEAST.Space, testqp::Matrix, trialqp::Matrix
 )
-    normals = normals(testqp)
+    trialnormals = normals(trialqp)
     rc, vals = potentials(testqp, testspace)
-    Btest = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
+    Btest = dropzeros(sparse(rc[:, 2], rc[:, 1], vals))
 
-    testspace == trialspace && return Btest, sparse(transpose(Btest)), normals
+    testspace == trialspace && return Btest, sparse(transpose(Btest)), trialnormals
     rc, vals = potentials(trialqp, trialspace)
-    Btrial = dropzeros(sparse(rc[:, 2], rc[:, 1], vals))
+    Btrial = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
 
-    return Btest, Btrial, normals
+    return Btest, Btrial, trialnormals
 end
diff --git a/ext/CFMMBEAST/operator/HH3Ddoublelayertransposed.jl b/ext/CFMMBEAST/operator/HH3Ddoublelayertransposed.jl
index 89c98ef..cf61f2a 100644
--- a/ext/CFMMBEAST/operator/HH3Ddoublelayertransposed.jl
+++ b/ext/CFMMBEAST/operator/HH3Ddoublelayertransposed.jl
@@ -5,11 +5,13 @@ function (fmmfunctor::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     testqp::Matrix,
     trialqp::Matrix,
     fmm;
-    ntasks=Threads.nthreads(),
+    scheduler,
 ) where {T<:BEAST.HH3DDoubleLayerTransposedFDBIO}
     Btest, Btrial, normals = HH3DDLTpotentialmatrix(testspace, trialspace, testqp, trialqp)
 
-    return CFMMMatrixHH3DDoubleLayerTransposed{scalartype(operator)}(
+    return CorrectionFactorMatrixMethod.CFMMMatrixHH3DDoubleLayerTransposed{
+        scalartype(operator)
+    }(
         operator, fmm, Btest, Btrial, normals
     )
 end
@@ -17,12 +19,12 @@ end
 function HH3DDLTpotentialmatrix(
     testspace::BEAST.Space, trialspace::BEAST.Space, testqp::Matrix, trialqp::Matrix
 )
-    normals = normals(testqp)
+    testnormals = normals(testqp)
     rc, vals = potentials(testqp, testspace)
-    Btest = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
-    testspace == trialspace && return Btest, sparse(transpose(Btest)), normals
+    Btest = dropzeros(sparse(rc[:, 2], rc[:, 1], vals))
+    testspace == trialspace && return Btest, sparse(transpose(Btest)), testnormals
 
-    rc_trial, vals_trial = potentials(trialqp, trialspace)
-    Btrial = dropzeros(sparse(rc_trial[:, 2], rc_trial[:, 1], vals_trial))
-    return Btest, Btrial, normals
+    rc, vals = potentials(trialqp, trialspace)
+    Btrial = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
+    return Btest, Btrial, testnormals
 end
diff --git a/ext/CFMMBEAST/operator/HH3Dhypersingular.jl b/ext/CFMMBEAST/operator/HH3Dhypersingular.jl
index 806356f..bc7f7e4 100644
--- a/ext/CFMMBEAST/operator/HH3Dhypersingular.jl
+++ b/ext/CFMMBEAST/operator/HH3Dhypersingular.jl
@@ -5,13 +5,13 @@ function (fmmfunctor::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     testqp::Matrix,
     trialqp::Matrix,
     fmm;
-    ntasks=Threads.nthreads(),
+    scheduler,
 ) where {T<:BEAST.HH3DHyperSingularFDBIO}
     testnormals, curlBtest, Btest, trialnormals, curlBtrial, Btrial = HH3DHSpotentialmatrix(
         testspace, trialspace, testqp, trialqp
     )
 
-    return FMMMatrixHS{scalartype(operator)}(
+    return CorrectionFactorMatrixMethod.CFMMMatrixHH3DHyperSingular{scalartype(operator)}(
         operator, fmm, testnormals, trialnormals, curlBtest, curlBtrial, Btest, Btrial
     )
 end
@@ -19,31 +19,31 @@ end
 function HH3DHSpotentialmatrix(
     testspace::BEAST.Space, trialspace::BEAST.Space, testqp::Matrix, trialqp::Matrix
 )
-    testnormals = getnormals(testqp)
+    testnormals = normals(testqp)
     rc_curl, vals_curl = curlpotentials(testqp, testspace)
     curlBtest = [
-        dropzeros(sparse(rc_curl[:, 1], rc_curl[:, 2], vals_curl[:, 1])),
-        dropzeros(sparse(rc_curl[:, 1], rc_curl[:, 2], vals_curl[:, 2])),
-        dropzeros(sparse(rc_curl[:, 1], rc_curl[:, 2], vals_curl[:, 3])),
+        dropzeros(sparse(rc_curl[:, 2], rc_curl[:, 1], vals_curl[:, 1])),
+        dropzeros(sparse(rc_curl[:, 2], rc_curl[:, 1], vals_curl[:, 2])),
+        dropzeros(sparse(rc_curl[:, 2], rc_curl[:, 1], vals_curl[:, 3])),
     ]
 
     rc, vals = potentials(testqp, testspace)
-    Btest = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
+    Btest = dropzeros(sparse(rc[:, 2], rc[:, 1], vals))
 
     testspace == trialspace && return testnormals,
-    curlBtest, Btest, testnormals, sparse.(transpose(curlBtest)),
+    curlBtest, Btest, testnormals, sparse.(transpose.(curlBtest)),
     sparse(transpose(Btest))
 
-    trialnormals = getnormals(trialqp)
+    trialnormals = normals(trialqp)
     rc_curl, vals_curl = curlpotentials(trialqp, trialspace)
     curlBtrial = [
-        dropzeros(sparse(rc_curl[:, 2], rc_curl[:, 1], vals_curl[:, 1])),
-        dropzeros(sparse(rc_curl[:, 2], rc_curl[:, 1], vals_curl[:, 2])),
-        dropzeros(sparse(rc_curl[:, 2], rc_curl[:, 1], vals_curl[:, 3])),
+        dropzeros(sparse(rc_curl[:, 1], rc_curl[:, 2], vals_curl[:, 1])),
+        dropzeros(sparse(rc_curl[:, 1], rc_curl[:, 2], vals_curl[:, 2])),
+        dropzeros(sparse(rc_curl[:, 1], rc_curl[:, 2], vals_curl[:, 3])),
     ]
 
     rc, vals = potentials(trialqp, trialspace)
-    Btrial = dropzeros(sparse(rc[:, 2], rc[:, 1], vals))
+    Btrial = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
 
     return testnormals, curlBtest, Btest, trialnormals, curlBtrial, Btrial
 end
diff --git a/ext/CFMMBEAST/operator/HH3Dsinglelayer.jl b/ext/CFMMBEAST/operator/HH3Dsinglelayer.jl
index 86cfbe7..4ebc700 100644
--- a/ext/CFMMBEAST/operator/HH3Dsinglelayer.jl
+++ b/ext/CFMMBEAST/operator/HH3Dsinglelayer.jl
@@ -5,7 +5,7 @@ function (fmmfunctor::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     testqp::Matrix,
     trialqp::Matrix,
     fmm;
-    ntasks=Threads.nthreads(),
+    scheduler,
 ) where {T<:BEAST.HH3DSingleLayerFDBIO}
     Btest, Btrial = HH3DSLpotentialmatrix(testspace, trialspace, testqp, trialqp)
 
@@ -18,11 +18,11 @@ function HH3DSLpotentialmatrix(
     testspace::BEAST.Space, trialspace::BEAST.Space, testqp::Matrix, trialqp::Matrix
 )
     rc, vals = potentials(testqp, testspace)
-    Btest = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
+    Btest = dropzeros(sparse(rc[:, 2], rc[:, 1], vals))
 
     testspace == trialspace && return Btest, sparse(transpose(Btest))
 
-    rc_test, vals_test = potentials(testqp, testspace)
-    Btrial = dropzeros(sparse(rc_test[:, 2], rc_test[:, 1], vals_test))
+    rc, vals = potentials(trialqp, trialspace)
+    Btrial = dropzeros(sparse(rc[:, 1], rc[:, 2], vals))
     return Btest, Btrial
 end
diff --git a/ext/CFMMBEAST/operator/MW3Ddoublelayer.jl b/ext/CFMMBEAST/operator/MW3Ddoublelayer.jl
index 09dad24..1ad441e 100644
--- a/ext/CFMMBEAST/operator/MW3Ddoublelayer.jl
+++ b/ext/CFMMBEAST/operator/MW3Ddoublelayer.jl
@@ -6,11 +6,13 @@ function (fmmfunctor::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     testqp::Matrix,
     trialqp::Matrix,
     fmm;
-    ntasks=Threads.nthreads(),
+    scheduler,
 ) where {T<:BEAST.MWDoubleLayer3D}
     Btest, Btrial = MW3DDLpotentialmatrix(testspace, trialspace, testqp, trialqp)
 
-    return CFMMMatrixMW3DDoubleLayer{scalartype(operator)}(operator, fmm, Btest, Btrial)
+    return CorrectionFactorMatrixMethod.CFMMMatrixMW3DDoubleLayer{scalartype(operator)}(
+        operator, fmm, Btest, Btrial
+    )
 end
 
 function MW3DDLpotentialmatrix(
@@ -18,17 +20,17 @@ function MW3DDLpotentialmatrix(
 )
     rc, vals = potentials(testqp, testspace)
     Btest = [
-        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 1])),
-        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 2])),
-        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 3])),
+        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 1])),
+        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 2])),
+        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 3])),
     ]
 
     testspace == trialspace && return Btest, sparse.(transpose.(Btest))
     rc, vals = potentials(trialqp, trialspace)
     Btrial = [
-        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 1])),
-        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 2])),
-        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 3])),
+        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 1])),
+        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 2])),
+        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 3])),
     ]
 
     return Btest, Btrial
diff --git a/ext/CFMMBEAST/operator/MW3Dsinglelayer.jl b/ext/CFMMBEAST/operator/MW3Dsinglelayer.jl
index b3c5a22..78289e5 100644
--- a/ext/CFMMBEAST/operator/MW3Dsinglelayer.jl
+++ b/ext/CFMMBEAST/operator/MW3Dsinglelayer.jl
@@ -6,13 +6,13 @@ function (::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     testqp::Matrix,
     trialqp::Matrix,
     fmm;
-    ntasks=Threads.nthreads(),
+    scheduler,
 ) where {T<:BEAST.MWSingleLayer3D}
     Btest, divBtest, Btrial, divBtrial = MW3DSLpotentialmatrix(
         testspace, trialspace, testqp, trialqp
     )
 
-    return CFMMMatrixMW3DSingleLayer{scalartype(operator)}(
+    return CorrectionFactorMatrixMethod.CFMMMatrixMW3DSingleLayer{scalartype(operator)}(
         operator, fmm, Btest, Btrial, divBtest, divBtrial
     )
 end
@@ -22,25 +22,25 @@ function MW3DSLpotentialmatrix(
 )
     rc, vals = potentials(testqp, testspace)
     Btest = [
-        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 1])),
-        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 2])),
-        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 3])),
+        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 1])),
+        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 2])),
+        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 3])),
     ]
 
     rcdiv, valsdiv = divpotentials(testqp, testspace)
-    divBtest = dropzeros(sparse(rcdiv[:, 1], rcdiv[:, 2], valsdiv))
+    divBtest = dropzeros(sparse(rcdiv[:, 2], rcdiv[:, 1], valsdiv))
 
     testspace == trialspace &&
         return Btest, divBtest, sparse.(transpose.(Btest)), sparse(transpose(divBtest))
 
     rc, vals = potentials(trialqp, trialspace)
     Btrial = [
-        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 1])),
-        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 2])),
-        dropzeros(sparse(rc[:, 2], rc[:, 1], vals[:, 3])),
+        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 1])),
+        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 2])),
+        dropzeros(sparse(rc[:, 1], rc[:, 2], vals[:, 3])),
     ]
     rcdiv, valsdiv = divpotentials(trialqp, trialspace)
-    divBtrial = dropzeros(sparse(rcdiv[:, 2], rcdiv[:, 1], valsdiv))
+    divBtrial = dropzeros(sparse(rcdiv[:, 1], rcdiv[:, 2], valsdiv))
 
     return Btest, divBtest, Btrial, divBtrial
 end
diff --git a/ext/CFMMBEAST/quadstrat.jl b/ext/CFMMBEAST/quadstrat.jl
index 896b5f3..61f9e40 100644
--- a/ext/CFMMBEAST/quadstrat.jl
+++ b/ext/CFMMBEAST/quadstrat.jl
@@ -22,7 +22,7 @@ function BEAST.quaddata(
     trial_elements,
     qs::SafeDoubleNumQStrat,
 )
-    test_quad_data  = BEAST.quadpoints(local_test_basis, test_elements, (qs.outer_rule,))
+    test_quad_data = BEAST.quadpoints(local_test_basis, test_elements, (qs.outer_rule,))
     trial_quad_data = BEAST.quadpoints(local_trial_basis, trial_elements, (qs.inner_rule,))
 
     return test_quad_data, trial_quad_data
diff --git a/ext/CFMMExaFMMt/CFMMExaFMMt.jl b/ext/CFMMExaFMMt/CFMMExaFMMt.jl
index 3ceda92..f9b5701 100644
--- a/ext/CFMMExaFMMt/CFMMExaFMMt.jl
+++ b/ext/CFMMExaFMMt/CFMMExaFMMt.jl
@@ -1,7 +1,6 @@
 module CFMMExaFMMt
 
 using CorrectionFactorMatrixMethod
-import CorrectionFactorMatrixMethod: ExaFMMtFunctor
 using BEAST
 using ExaFMMt
 
@@ -11,4 +10,8 @@ function CorrectionFactorMatrixMethod.setup(spoints, tpoints, options::ExaFMMt.F
     return ExaFMMt.setup(spoints, tpoints, options)
 end
 
+function CorrectionFactorMatrixMethod.fmmresult(fmm::ExaFMMt.ExaFMM, x)
+    return ExaFMMt.evaluate(fmm, x, fmm.fmmoptions)
+end
+
 end
diff --git a/ext/CFMMExaFMMt/fmmoptions.jl b/ext/CFMMExaFMMt/fmmoptions.jl
index 8f4e493..93209ee 100644
--- a/ext/CFMMExaFMMt/fmmoptions.jl
+++ b/ext/CFMMExaFMMt/fmmoptions.jl
@@ -1,12 +1,19 @@
-exafmmoptions(
+function exafmmoptions(
     gamma::T, fmm::CorrectionFactorMatrixMethod.ExaFMMtFunctor
-) where {T<:Val{0}} = LaplaceFMMOptions(; p=fmm.p, ncrit=fmm.ncrit)
+) where {T<:Val{0}}
+    return LaplaceFMMOptions(; p=fmm.p, ncrit=fmm.ncrit)
+end
 #TODO: Write unit tests for the ModifiedHelmholtzFMMOptions
-exafmmoptions(gamma::T, fmm::CorrectionFactorMatrixMethod.ExaFMMtFunctor) where {T<:Real} =
-    ModifiedHelmholtzFMMOptions(gamma; p=fmm.p, ncrit=fmm.ncrit)
-exafmmoptions(
+function exafmmoptions(
+    gamma::T, fmm::CorrectionFactorMatrixMethod.ExaFMMtFunctor
+) where {T<:Real}
+    return ModifiedHelmholtzFMMOptions(gamma; p=fmm.p, ncrit=fmm.ncrit)
+end
+function exafmmoptions(
     gamma::T, fmm::CorrectionFactorMatrixMethod.ExaFMMtFunctor
-) where {T<:Complex} = HelmholtzFMMOptions(-gamma / im; p=fmm.p, ncrit=fmm.ncrit)
+) where {T<:Complex}
+    return HelmholtzFMMOptions(-gamma / im; p=fmm.p, ncrit=fmm.ncrit)
+end
 
 function (functor::CorrectionFactorMatrixMethod.ExaFMMtFunctor)(
     operator::BEAST.IntegralOperator
diff --git a/src/CFMM/assemble.jl b/src/CFMM/assemble.jl
new file mode 100644
index 0000000..3230458
--- /dev/null
+++ b/src/CFMM/assemble.jl
@@ -0,0 +1,66 @@
+function defaultminhalfsize(testspace::AbstractVector, trialspace::AbstractVector)
+    _, testhalfsize = H2Trees.boundingbox(testspace)
+    _, trialhalfsize = H2Trees.boundingbox(trialspace)
+    return max(testhalfsize, trialhalfsize) / 2^10
+end
+
+defaultminvalues(::FMMFunctor) = 50
+defaultminvalues(functor::ExaFMMtFunctor) = functor.ncrit
+
+"""
+    CFMM
+
+High-level assembly interface for correction-factor FMM operators.
+"""
+module CFMM
+
+using H2Trees
+using ..CorrectionFactorMatrixMethod:
+    FMMFunctor, PetrovGalerkinCFMM, defaultminhalfsize, defaultminvalues
+
+"""
+    assemble(operator, testspace, trialspace; kwargs...)
+    assemble(operator, space; kwargs...)
+
+Assemble a matrix-free correction-factor FMM operator.
+
+By default, an `H2Trees.TwoNTree` is constructed with a scale-dependent minimum
+box size and a target leaf occupancy equal to the FMM `ncrit` setting. Pass
+`tree` to use an existing tree, or override `minhalfsize`, `testminvalues`, and
+`trialminvalues`.
+
+All remaining keywords are forwarded to [`PetrovGalerkinCFMM`](@ref).
+"""
+function assemble(
+    operator,
+    testspace,
+    trialspace;
+    fmmfunctor=FMMFunctor(),
+    tree=nothing,
+    minhalfsize=nothing,
+    testminvalues=defaultminvalues(fmmfunctor),
+    trialminvalues=defaultminvalues(fmmfunctor),
+    treeoptions=(;),
+    kwargs...,
+)
+    if isnothing(tree)
+        isnothing(minhalfsize) && (minhalfsize = defaultminhalfsize(testspace, trialspace))
+        tree = H2Trees.TwoNTree(
+            testspace,
+            trialspace,
+            minhalfsize;
+            testminvalues=testminvalues,
+            trialminvalues=trialminvalues,
+            treeoptions...,
+        )
+    end
+    return PetrovGalerkinCFMM(
+        operator, testspace, trialspace, tree; fmmfunctor=fmmfunctor, kwargs...
+    )
+end
+
+function assemble(operator, space; kwargs...)
+    return assemble(operator, space, space; kwargs...)
+end
+
+end
diff --git a/src/CFMM/galerkinfmm.jl b/src/CFMM/galerkinfmm.jl
index e69de29..8b13789 100644
--- a/src/CFMM/galerkinfmm.jl
+++ b/src/CFMM/galerkinfmm.jl
@@ -0,0 +1 @@
+
diff --git a/src/CFMM/petrovgalerkincfmm.jl b/src/CFMM/petrovgalerkincfmm.jl
index 32f4ff4..6327120 100644
--- a/src/CFMM/petrovgalerkincfmm.jl
+++ b/src/CFMM/petrovgalerkincfmm.jl
@@ -1,11 +1,14 @@
-struct PetrovGalerkinCFMM{K,CorrectedNearsType,FMMType} <: LinearMaps.LinearMap{K}
+struct PetrovGalerkinCFMM{K,CorrectedNearsType,FMMType,SchedulerType} <:
+       LinearMaps.LinearMap{K}
     correctednears::CorrectedNearsType
     fmm::FMMType
     dim::Tuple{Int,Int}
-    ntasks::Int
+    scheduler::SchedulerType
 
-    function PetrovGalerkinCFMM{K}(correctednears, fmm, dim, ntasks) where {K}
-        return new{K,typeof(correctednears),typeof(fmm)}(correctednears, fmm, dim, ntasks)
+    function PetrovGalerkinCFMM{K}(correctednears, fmm, dim, scheduler) where {K}
+        return new{K,typeof(correctednears),typeof(fmm),typeof(scheduler)}(
+            correctednears, fmm, dim, scheduler
+        )
     end
 end
 
@@ -15,6 +18,17 @@ end
 
 function scalartype(operator) end
 
+"""
+    PetrovGalerkinCFMM(operator, testspace, trialspace, tree; kwargs...)
+
+Construct a matrix-free Petrov-Galerkin operator whose far interactions are
+evaluated by a fast multipole method and whose near interactions are corrected
+by direct quadrature.
+
+The concrete operator and space methods are supplied by package extensions.
+Loading `BEAST`, `CompScienceMeshes`, and `ExaFMMt` enables the supported
+boundary-element implementation.
+"""
 function PetrovGalerkinCFMM(
     operator,
     testspace,
@@ -23,65 +37,57 @@ function PetrovGalerkinCFMM(
     fmmfunctor=FMMFunctor(),
     nearquadstrat=defaultnearquadstrat(operator, testspace, trialspace),
     farquadstrat=defaultfarquadstrat(operator, testspace, trialspace),
-    ntasks=Threads.nthreads(),
+    scheduler=DynamicScheduler(),
     isnear=H2Trees.isnear,
     computetransposeadjoint=false,
 )
-    correctedkernelmatrix = AbstractCorrectedKernelMatrix(
+    correctednears = assemblecorrectednears(
         operator,
         testspace,
-        trialspace;
+        trialspace,
+        tree;
         nearquadstrat=nearquadstrat,
         farquadstrat=farquadstrat,
-    )
-    values, nearvalues = nearinteractions(tree; isnear=isnear)
-    blocks = tmap(values, nearvalues) do v, nv
-        blk = zeros(scalartype(operator), length(v), length(nv))
-        correctedkernelmatrix(v, nv, blk)
-        return blk
-    end
-    correctednears = BlockSparseMatrix(
-        blocks, values, nearvalues, length(testspace), length(trialspace)
+        scheduler=scheduler,
+        isnear=isnear,
     )
 
     fmm = fmmfunctor(
         operator,
         testspace,
         trialspace;
-        ntasks=ntasks,
+        scheduler=scheduler,
         farquadstrat=farquadstrat,
         computetransposeadjoint=computetransposeadjoint,
     )
 
     return PetrovGalerkinCFMM{scalartype(operator)}(
-        correctednears, fmm, (length(testspace), length(trialspace)), ntasks
+        correctednears, fmm, (length(testspace), length(trialspace)), scheduler
     )
 end
 
-function LinearMaps.mul!(
-    y::AbstractVector{K}, A::PetrovGalerkinCFMM{K}, x::AbstractVector{K}
-) where {K}
+function LinearMaps.mul!(y::AbstractVector, A::PetrovGalerkinCFMM, x::AbstractVector)
     LinearMaps.mul!(y, A.fmm, x)
-    y += A.correctednears * x
+    mul!(y, A.correctednears, x, true, true)
     return y
 end
 
 function LinearMaps.mul!(
-    y::AbstractVector{K},
-    A::LinearMaps.TransposeMap{<:Any,<:PetrovGalerkinCFMM{K}},
-    x::AbstractVector{K},
-) where {K}
+    y::AbstractVector,
+    A::LinearMaps.TransposeMap{<:Any,<:PetrovGalerkinCFMM},
+    x::AbstractVector,
+)
     LinearMaps.mul!(y, transpose(A.lmap.fmm), x)
-    y += transpose(A.lmap.correctednears) * x
+    mul!(y, transpose(A.lmap.correctednears), x, true, true)
     return y
 end
 
 function LinearMaps.mul!(
-    y::AbstractVector{K},
-    A::LinearMaps.AdjointMap{<:Any,<:PetrovGalerkinCFMM{K}},
-    x::AbstractVector{K},
-) where {K}
+    y::AbstractVector,
+    A::LinearMaps.AdjointMap{<:Any,<:PetrovGalerkinCFMM},
+    x::AbstractVector,
+)
     LinearMaps.mul!(y, adjoint(A.lmap.fmm), x)
-    y += adjoint(A.lmap.correctednears) * x
+    mul!(y, adjoint(A.lmap.correctednears), x, true, true)
     return y
 end
diff --git a/src/CorrectionFactorMatrixMethod.jl b/src/CorrectionFactorMatrixMethod.jl
index b829a53..f543f8f 100644
--- a/src/CorrectionFactorMatrixMethod.jl
+++ b/src/CorrectionFactorMatrixMethod.jl
@@ -1,3 +1,25 @@
+"""
+    CorrectionFactorMatrixMethod
+
+Matrix-free boundary-element operators that combine a fast multipole
+approximation of the far interactions with directly assembled near-field
+corrections.
+
+The far field is evaluated by a fast multipole method, while the near
+interactions are assembled with the boundary-element quadrature and corrected
+for the part already represented by the FMM:
+
+    A x ≈ A_fmm x + (A_near - A_fmm,near) x
+
+The generic correction and linear-map machinery lives in the core package;
+support for concrete operators and spaces is provided through package
+extensions for [BEAST](https://github.com/krcools/BEAST.jl),
+[CompScienceMeshes](https://github.com/krcools/CompScienceMeshes.jl), and
+[ExaFMMt](https://github.com/JoshuaTetzner/ExaFMMt.jl).
+
+The main entry point is [`CFMM.assemble`](@ref); see also
+[`PetrovGalerkinCFMM`](@ref) and [`FMMFunctor`](@ref).
+"""
 module CorrectionFactorMatrixMethod
 
 using BlockSparseMatrices
@@ -42,7 +64,8 @@ include("nearinteractions.jl")
 include("operators/abstractoperators.jl")
 
 include("CFMM/petrovgalerkincfmm.jl")
+include("CFMM/assemble.jl")
 
-export PetrovGalerkinCFMM
+export CFMM, PetrovGalerkinCFMM
 
 end
diff --git a/src/kernelmatrix/abstractcorrectedkernelmatrix.jl b/src/kernelmatrix/abstractcorrectedkernelmatrix.jl
index 4e3f9cc..076c32a 100644
--- a/src/kernelmatrix/abstractcorrectedkernelmatrix.jl
+++ b/src/kernelmatrix/abstractcorrectedkernelmatrix.jl
@@ -2,5 +2,5 @@ abstract type AbstractCorrectedKernelMatrix{T} end
 
 function AbstractCorrectedKernelMatrix(operator, testspace, trialspace; args...) end
 
-function (::AbstractCorrectedKernelMatrix)(tdata, sdata, matrixblock) end
+function (::AbstractCorrectedKernelMatrix)(matrixblock, tdata, sdata) end
 Base.eltype(::AbstractCorrectedKernelMatrix{T}) where {T} = T
diff --git a/src/nearinteractions.jl b/src/nearinteractions.jl
index ad6b95e..539a0f4 100644
--- a/src/nearinteractions.jl
+++ b/src/nearinteractions.jl
@@ -1,47 +1,33 @@
-import H2Trees: isleaf, testtree, trialtree, root, children
-
-function nears!(
-    tree,
-    values::Vector{V},
-    nearvalues::Vector{V},
-    tnode::Int,
-    snodes::V;
+function assemblecorrectednears(
+    operator,
+    testspace,
+    trialspace,
+    tree;
+    nearquadstrat=defaultnearquadstrat(operator, testspace, trialspace),
+    farquadstrat=defaultfarquadstrat(operator, testspace, trialspace),
+    scheduler=DynamicScheduler(),
     isnear=H2Trees.isnear,
-) where {V<:Vector{Int}}
-    nearnodes = Int[]
-    childnearnodes = Int[]
-    for snode in snodes
-        if isnear(testtree(tree), trialtree(tree), tnode, snode)
-            if isleaf(testtree(tree), tnode) || isleaf(trialtree(tree), snode)
-                push!(nearnodes, snode)
-            else
-                append!(childnearnodes, collect(children(trialtree(tree), snode)))
-            end
-        end
-    end
-    if nearnodes != []
-        push!(nearvalues, H2Trees.values(trialtree(tree), nearnodes))
-        push!(values, H2Trees.values(testtree(tree), tnode))
-    end
-    if childnearnodes != []
-        for child in children(testtree(tree), tnode)
-            nears!(tree, values, nearvalues, child, childnearnodes; isnear=isnear)
-        end
+)
+    matrix = AbstractCorrectedKernelMatrix(
+        operator,
+        testspace,
+        trialspace;
+        nearquadstrat=nearquadstrat,
+        farquadstrat=farquadstrat,
+    )
+    values, nearvalues = H2Trees.nearinteractions(tree; isnear=isnear)
+    blocks = tmap(values, nearvalues; scheduler=scheduler) do v, nv
+        block = zeros(scalartype(operator), length(v), length(nv))
+        matrix(block, v, nv)
+        return block
     end
-end
-
-function nearinteractions(tree::H2Trees.BlockTree; isnear=H2Trees.isnear)
-    !isnear(testtree(tree), trialtree(tree), root(testtree(tree)), root(trialtree(tree))) &&
-        return Vector{Int}(), Vector{Int}[]
-    values = Vector{Int}[]
-    nearvalues = Vector{Int}[]
-    nears!(
-        tree,
+    scheduler = isempty(blocks) ? SerialScheduler() : scheduler
+    return BlockSparseMatrix(
+        blocks,
         values,
         nearvalues,
-        root(testtree(tree)),
-        [root(trialtree(tree))];
-        isnear=isnear,
+        length(testspace),
+        length(trialspace);
+        scheduler=scheduler,
     )
-    return values, nearvalues
 end
diff --git a/src/operators/HH3Ddoublelayer.jl b/src/operators/HH3Ddoublelayer.jl
index 9e02bd2..a3bf28d 100644
--- a/src/operators/HH3Ddoublelayer.jl
+++ b/src/operators/HH3Ddoublelayer.jl
@@ -43,9 +43,9 @@ end
         xfmm = x
     end
 
-    fmm_res1 = A.Btest * (A.fmm * (A.normals[:, 1] .* (A.Btrial * xfmm)))[:, 2]
-    fmm_res2 = A.Btest * (A.fmm * (A.normals[:, 2] .* (A.Btrial * xfmm)))[:, 3]
-    fmm_res3 = A.Btest * (A.fmm * (A.normals[:, 3] .* (A.Btrial * xfmm)))[:, 4]
+    fmm_res1 = A.Btest * fmmresult(A.fmm, A.normals[:, 1] .* (A.Btrial * xfmm))[:, 2]
+    fmm_res2 = A.Btest * fmmresult(A.fmm, A.normals[:, 2] .* (A.Btrial * xfmm))[:, 3]
+    fmm_res3 = A.Btest * fmmresult(A.fmm, A.normals[:, 3] .* (A.Btrial * xfmm))[:, 4]
     fmm_res = -(fmm_res1 + fmm_res2 + fmm_res3)
     y .= A.operator.alpha .* fmm_res
 
@@ -75,7 +75,7 @@ end
         xfmm = x
     end
 
-    xx = -1 .* (transpose(A.fmm) * (transpose(A.Btest) * xfmm))
+    xx = -1 .* fmmresult(transpose(A.fmm), transpose(A.Btest) * xfmm)
     test = transpose(A.Btrial)
 
     fmm_res1 = test * (A.normals[:, 1] .* xx[:, 2])
diff --git a/src/operators/HH3Ddoublelayertransposed.jl b/src/operators/HH3Ddoublelayertransposed.jl
index 92f5f31..24e7d5a 100644
--- a/src/operators/HH3Ddoublelayertransposed.jl
+++ b/src/operators/HH3Ddoublelayertransposed.jl
@@ -44,7 +44,7 @@ end
         xfmm = x
     end
 
-    res = A.fmm.A * (A.Btrial * xfmm)
+    res = fmmresult(A.fmm, A.Btrial * xfmm)
     fmm_res1 = A.normals[:, 1] .* (res)[:, 2]
     fmm_res2 = A.normals[:, 2] .* (res)[:, 3]
     fmm_res3 = A.normals[:, 3] .* (res)[:, 4]
@@ -78,9 +78,9 @@ end
     end
 
     xx = transpose(A.Btest) * xfmm
-    fmm_res1 = (transpose(A.fmm) * (A.normals[:, 1] .* xx))[:, 2]
-    fmm_res2 = (transpose(A.fmm) * (A.normals[:, 2] .* xx))[:, 3]
-    fmm_res3 = (transpose(A.fmm) * (A.normals[:, 3] .* xx))[:, 4]
+    fmm_res1 = fmmresult(transpose(A.fmm), A.normals[:, 1] .* xx)[:, 2]
+    fmm_res2 = fmmresult(transpose(A.fmm), A.normals[:, 2] .* xx)[:, 3]
+    fmm_res3 = fmmresult(transpose(A.fmm), A.normals[:, 3] .* xx)[:, 4]
     fmm_res = -(fmm_res1 + fmm_res2 + fmm_res3)
 
     y .= A.op.alpha .* (transpose(A.Btrial) * fmm_res)
diff --git a/src/operators/HH3Dhypersingular.jl b/src/operators/HH3Dhypersingular.jl
index 5ca311a..0d5eaaf 100644
--- a/src/operators/HH3Dhypersingular.jl
+++ b/src/operators/HH3Dhypersingular.jl
@@ -19,7 +19,7 @@ struct CFMMMatrixHH3DHyperSingular{K,OperatorType,FMMType,NormalsType,SparseMatr
             typeof(testnormals),
             typeof(Btest),
         }(
-            op, fmm, testnormals, trialnormals, curlBtest, curlBtrial, Btest, Btrial
+            operator, fmm, testnormals, trialnormals, curlBtest, curlBtrial, Btest, Btrial
         )
     end
 end
@@ -50,7 +50,7 @@ end
     else
         xfmm = x
     end
-    if A.op.alpha != 0.0
+    if A.operator.alpha != 0.0
         fmm_res1 =
             A.testnormals[:, 1] .*
             (A.fmm * (A.trialnormals[:, 1] .* (A.Btrial * xfmm)))[:, 1]
diff --git a/src/operators/HH3Dsinglelayer.jl b/src/operators/HH3Dsinglelayer.jl
index 10744f7..42ed1c9 100644
--- a/src/operators/HH3Dsinglelayer.jl
+++ b/src/operators/HH3Dsinglelayer.jl
@@ -58,7 +58,7 @@ end
     eltype(x) != eltype(A.fmm) ? (xfmm = eltype(A.fmm).(x)) : (xfmm = x)
 
     y .=
-        A.op.alpha .*
+        alpha(A.operator) .*
         (transpose(A.Btrial) * (transpose(A.fmm) * (transpose(A.Btest) * xfmm))[:, 1])
 
     return y
diff --git a/src/operators/MW3Ddoublelayer.jl b/src/operators/MW3Ddoublelayer.jl
index 6117688..7b3c09b 100644
--- a/src/operators/MW3Ddoublelayer.jl
+++ b/src/operators/MW3Ddoublelayer.jl
@@ -7,7 +7,7 @@ struct CFMMMatrixMW3DDoubleLayer{K,OperatorType,FMMType,SparseMatrixType} <:
     Btrial::Vector{SparseMatrixType}
 
     function CFMMMatrixMW3DDoubleLayer{K}(operator, fmm, Btest, Btrial) where {K}
-        return new{scalartype(operator),typeof(operator),typeof(fmm),typeof(Btest)}(
+        return new{scalartype(operator),typeof(operator),typeof(fmm),eltype(Btest)}(
             operator, fmm, Btest, Btrial
         )
     end
@@ -31,9 +31,9 @@ end
     end
     fill!(y, zero(K))
 
-    res1 = (A.fmm * (A.Btrial[1] * x))[:, 2:4]
-    res2 = (A.fmm * (A.Btrial[2] * x))[:, 2:4]
-    res3 = (A.fmm * (A.Btrial[3] * x))[:, 2:4]
+    res1 = fmmresult(A.fmm, A.Btrial[1] * x)[:, 2:4]
+    res2 = fmmresult(A.fmm, A.Btrial[2] * x)[:, 2:4]
+    res3 = fmmresult(A.fmm, A.Btrial[3] * x)[:, 2:4]
 
     y1 = A.Btest[1] * (res3[:, 2] - res2[:, 3])
     y2 = A.Btest[2] * (res1[:, 3] - res3[:, 1])
@@ -57,9 +57,9 @@ end
     end
     fill!(y, zero(K))
 
-    res1 = (transpose(A.fmm) * (transpose(A.Btest[1]) * x))[:, 2:4]
-    res2 = (transpose(A.fmm) * (transpose(A.Btest[2]) * x))[:, 2:4]
-    res3 = (transpose(A.fmm) * (transpose(A.Btest[3]) * x))[:, 2:4]
+    res1 = fmmresult(transpose(A.fmm), transpose(A.Btest[1]) * x)[:, 2:4]
+    res2 = fmmresult(transpose(A.fmm), transpose(A.Btest[2]) * x)[:, 2:4]
+    res3 = fmmresult(transpose(A.fmm), transpose(A.Btest[3]) * x)[:, 2:4]
 
     y1 = transpose(A.Btrial[1]) * (res3[:, 2] - res2[:, 3])
     y2 = transpose(A.Btrial[2]) * (res1[:, 3] - res3[:, 1])
diff --git a/src/operators/MW3Dsinglelayer.jl b/src/operators/MW3Dsinglelayer.jl
index e1a0cca..8f7f3f7 100644
--- a/src/operators/MW3Dsinglelayer.jl
+++ b/src/operators/MW3Dsinglelayer.jl
@@ -11,7 +11,7 @@ struct CFMMMatrixMW3DSingleLayer{K,OperatorType,FMMType,SparseMatrixType} <:
     function CFMMMatrixMW3DSingleLayer{K}(
         operator, fmm, Btest, Btrial, divBtest, divBtrial
     ) where {K}
-        return new{scalartype(operator),typeof(operator),typeof(fmm),typeof(Btest)}(
+        return new{scalartype(operator),typeof(operator),typeof(fmm),eltype(Btest)}(
             operator, fmm, Btest, Btrial, divBtest, divBtrial
         )
     end
diff --git a/src/operators/abstractoperators.jl b/src/operators/abstractoperators.jl
index 0b27242..945fa83 100644
--- a/src/operators/abstractoperators.jl
+++ b/src/operators/abstractoperators.jl
@@ -16,6 +16,20 @@ struct ExaFMMtFunctor <: FMMFunctor
     ncrit::Int
 end
 
+"""
+    FMMFunctor(; p=8, ncrit=50)
+
+Configure the fast multipole method used for the far interactions.
+
+`p` is the expansion order and `ncrit` the critical leaf size, i.e. the maximum
+number of quadrature points per box. A larger `p` increases accuracy and cost,
+while `ncrit` trades tree depth against the size of the direct near-field
+blocks.
+
+The returned functor is passed as the `fmmfunctor` keyword to
+[`CFMM.assemble`](@ref) and [`PetrovGalerkinCFMM`](@ref); its default targets
+the ExaFMMt backend.
+"""
 FMMFunctor(; p=8, ncrit=50) = ExaFMMtFunctor(p, ncrit)
 
 function (functor::ExaFMMtFunctor)(operator)
@@ -32,16 +46,22 @@ struct FMM{K,FMMType,TransposeFMMType} <: LinearMaps.LinearMap{K}
     end
 end
 Base.eltype(fmm::FMM) = eltype(fmm.A)
-Base.size(fmm::FMM, dim=nothing) = size(fmm.A, dim)
+Base.size(fmm::FMM) = size(fmm.A)
+Base.size(fmm::FMM, dim::Integer) = size(fmm.A, dim)
 
-LinearMaps._unsafe_mul!(y::AbstractVecOrMat, fmm::FMM, x::AbstractVector) =
-    mul!(y, fmm.A, x)
-LinearMaps._unsafe_mul!(
+function LinearMaps._unsafe_mul!(y::AbstractVecOrMat, fmm::FMM, x::AbstractVector)
+    return mul!(y, fmm.A, x)
+end
+function LinearMaps._unsafe_mul!(
     y::AbstractVecOrMat, fmm::LinearMaps.TransposeMap{<:Any,<:FMM}, x::AbstractVector
-) = mul!(y, fmm.lmap.Aᵀ, x)
-LinearMaps._unsafe_mul!(
+)
+    return mul!(y, fmm.lmap.Aᵀ, x)
+end
+function LinearMaps._unsafe_mul!(
     y::AbstractVecOrMat, fmm::LinearMaps.AdjointMap{<:Any,<:FMM}, x::AbstractVector
-) = error("Adjoint multiplication not implemented for FMM")
+)
+    return error("Adjoint multiplication not implemented for FMM")
+end
 struct SymmetricFMM{K,FMMType} <: LinearMaps.LinearMap{K}
     A::FMMType
 
@@ -50,18 +70,24 @@ struct SymmetricFMM{K,FMMType} <: LinearMaps.LinearMap{K}
     end
 end
 Base.eltype(fmm::SymmetricFMM) = eltype(fmm.A)
-Base.size(fmm::SymmetricFMM, dim=nothing) = size(fmm.A, dim)
+Base.size(fmm::SymmetricFMM) = size(fmm.A)
+Base.size(fmm::SymmetricFMM, dim::Integer) = size(fmm.A, dim)
 
-LinearMaps._unsafe_mul!(y::AbstractVecOrMat, fmm::SymmetricFMM, x::AbstractVector) =
-    mul!(y, fmm.A, x)
-LinearMaps._unsafe_mul!(
+function LinearMaps._unsafe_mul!(y::AbstractVecOrMat, fmm::SymmetricFMM, x::AbstractVector)
+    return mul!(y, fmm.A, x)
+end
+function LinearMaps._unsafe_mul!(
     y::AbstractVecOrMat,
     fmm::LinearMaps.TransposeMap{<:Any,<:SymmetricFMM},
     x::AbstractVector,
-) = mul!(y, fmm.lmap.A, x)
-LinearMaps._unsafe_mul!(
+)
+    return mul!(y, fmm.lmap.A, x)
+end
+function LinearMaps._unsafe_mul!(
     y::AbstractVecOrMat, fmm::LinearMaps.AdjointMap{<:Any,<:SymmetricFMM}, x::AbstractVector
-) = error("Adjoint multiplication not implemented for SymmetricFMM")
+)
+    return error("Adjoint multiplication not implemented for SymmetricFMM")
+end
 
 function FMM(A::T) where {T}
     return SymmetricFMM{eltype(A)}(A)
@@ -102,6 +128,13 @@ function setup(spoints, tpoints, options)
     return error("This function has to be implemented for ", typeof(options))
 end
 
+function fmmresult end
+
+fmmresult(fmm::FMM, x) = fmmresult(fmm.A, x)
+fmmresult(fmm::LinearMaps.TransposeMap{<:Any,<:FMM}, x) = fmmresult(fmm.lmap.Aᵀ, x)
+fmmresult(fmm::SymmetricFMM, x) = fmmresult(fmm.A, x)
+fmmresult(fmm::LinearMaps.TransposeMap{<:Any,<:SymmetricFMM}, x) = fmmresult(fmm.lmap.A, x)
+
 function setup_fmm(
     spoints::Matrix{F}, tpoints::Matrix{F}, options; computetransposeadjoint=false
 ) where {F<:Real}
@@ -116,17 +149,18 @@ function (fmmfunctor::FMMFunctor)(
     trialspace;
     farquadstrat=defaultfarquadstrat(operator, testspace, trialspace),
     computetransposeadjoint=false,
-    ntasks=Threads.nthreads(),
+    scheduler=DynamicScheduler(),
 )
     testsrcs, testqp = sources(testspace, farquadstrat.outer_rule)
     trialsrcs, trialqp = sources(trialspace, farquadstrat.inner_rule)
+    computetransposeadjoint |= testsrcs != trialsrcs
 
     fmm = setup_fmm(
-        testsrcs,
         trialsrcs,
+        testsrcs,
         fmmfunctor(operator);
         computetransposeadjoint=computetransposeadjoint,
     )
 
-    return fmmfunctor(operator, testspace, trialspace, testqp, trialqp, fmm; ntasks=ntasks)
+    return fmmfunctor(operator, testspace, trialspace, testqp, trialqp, fmm; scheduler)
 end
diff --git a/test/CFMMBEAST/test_CFMMBEAST.jl b/test/CFMMBEAST/test_CFMMBEAST.jl
new file mode 100644
index 0000000..6fd2fe1
--- /dev/null
+++ b/test/CFMMBEAST/test_CFMMBEAST.jl
@@ -0,0 +1,88 @@
+@testset "CFMMBEAST extension" begin
+    @test Base.get_extension(CFM, :CFMMBEAST) !== nothing
+
+    mesh = meshrectangle(1.0, 1.0, 0.5)
+    scalarspace = lagrangec0d1(mesh)
+    vectorspace = raviartthomas(mesh)
+    operator = Helmholtz3D.singlelayer(; wavenumber=0.5)
+    nearquadstrat = CFM.defaultnearquadstrat(operator, scalarspace, scalarspace)
+    farquadstrat = CFM.SafeDoubleNumQStrat(3, 3)
+
+    @test CFM.scalartype(operator) == BEAST.scalartype(operator)
+    @test CFM.alpha(operator) == operator.alpha
+    @test CFM.beta(Helmholtz3D.hypersingular(; wavenumber=0.5)) ==
+        Helmholtz3D.hypersingular(; wavenumber=0.5).beta
+
+    matrix = CFM.AbstractCorrectedKernelMatrix(
+        operator, scalarspace, scalarspace; nearquadstrat=nearquadstrat
+    )
+    block = zeros(ComplexF64, size(matrix))
+    matrix(block, collect(axes(block, 1)), collect(axes(block, 2)))
+    near = assemble(operator, scalarspace, scalarspace; quadstrat=nearquadstrat)
+    far = assemble(operator, scalarspace, scalarspace; quadstrat=farquadstrat)
+    @test block ≈ near - far
+
+    for space in (scalarspace, vectorspace)
+        points, quadrature = CFM.sources(space, 2)
+        indices, values = CFM.potentials(quadrature, space)
+        @test size(points, 2) == 3
+        @test size(indices, 2) == 2
+        @test size(values, 1) == size(indices, 1)
+    end
+
+    _, scalarquadrature = CFM.sources(scalarspace, 2)
+    curlindices, curlvalues = CFM.curlpotentials(scalarquadrature, scalarspace)
+    normalvalues = CFM.normals(scalarquadrature)
+    @test size(curlvalues) == (size(curlindices, 1), 3)
+    @test size(normalvalues, 2) == 3
+
+    _, vectorquadrature = CFM.sources(vectorspace, 2)
+    divindices, divvalues = CFM.divpotentials(vectorquadrature, vectorspace)
+    @test length(divvalues) == size(divindices, 1)
+
+    for (operator, space, matrixtype) in (
+        (
+            Helmholtz3D.singlelayer(; wavenumber=0.5),
+            scalarspace,
+            CFM.CFMMMatrixHH3DSingleLayer,
+        ),
+        (
+            Helmholtz3D.doublelayer(; wavenumber=0.5),
+            scalarspace,
+            CFM.CFMMMatrixHH3DDoubleLayer,
+        ),
+        (
+            Helmholtz3D.doublelayer_transposed(; wavenumber=0.5),
+            scalarspace,
+            CFM.CFMMMatrixHH3DDoubleLayerTransposed,
+        ),
+        (
+            Helmholtz3D.hypersingular(; wavenumber=0.5),
+            scalarspace,
+            CFM.CFMMMatrixHH3DHyperSingular,
+        ),
+        (
+            Maxwell3D.singlelayer(; wavenumber=0.5),
+            vectorspace,
+            CFM.CFMMMatrixMW3DSingleLayer,
+        ),
+        (
+            Maxwell3D.doublelayer(; wavenumber=0.5),
+            vectorspace,
+            CFM.CFMMMatrixMW3DDoubleLayer,
+        ),
+    )
+        _, quadrature = CFM.sources(space, 2)
+        matrix = CFM.FMMFunctor()(
+            operator,
+            space,
+            space,
+            quadrature,
+            quadrature,
+            nothing;
+            scheduler=CFM.SerialScheduler(),
+        )
+        @test matrix isa matrixtype
+        @test size(matrix) == (length(space), length(space))
+    end
+end
diff --git a/test/CFMMExaFMMt/test_CFMMExaFMMt.jl b/test/CFMMExaFMMt/test_CFMMExaFMMt.jl
new file mode 100644
index 0000000..a82abd9
--- /dev/null
+++ b/test/CFMMExaFMMt/test_CFMMExaFMMt.jl
@@ -0,0 +1,22 @@
+@testset "CFMMExaFMMt extension" begin
+    @test Base.get_extension(CFM, :CFMMExaFMMt) !== nothing
+
+    functor = CFM.FMMFunctor(; p=8, ncrit=32)
+    operator = Helmholtz3D.singlelayer(; wavenumber=0.5)
+    options = functor(operator)
+    @test options isa ExaFMMt.FMMOptions
+    @test (options.p, options.ncrit) == (8, 32)
+
+    sources = [0.0 0.0 0.0; 1.0 0.0 0.0]
+    targets = [0.0 1.0 0.0; 1.0 1.0 0.0]
+    fmm = CFM.setup(sources, targets, options)
+    strengths = ones(ComplexF64, size(sources, 1))
+    result = CFM.fmmresult(fmm, strengths)
+    wrapped = CFM.setup_fmm(sources, targets, options; computetransposeadjoint=true)
+
+    @test fmm isa ExaFMMt.ExaFMM
+    @test result == ExaFMMt.evaluate(fmm, strengths, fmm.fmmoptions)
+    @test size(result) == (size(targets, 1), 4)
+    @test wrapped isa CFM.FMM
+    @test CFM.fmmresult(wrapped, strengths) == result
+end
diff --git a/test/blockassembler.jl b/test/blockassembler.jl
deleted file mode 100644
index 4794cce..0000000
--- a/test/blockassembler.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-using BEAST
-using CompScienceMeshes
-using H2Trees
-using ExaFMMt
-using CorrectionFactorMatrixMethod
-
-Γ = meshsphere(1.0, 0.1)
-
-op = Helmholtz3D.singlelayer()
-space = lagrangecxd0(Γ)
-tree = H2Trees.TwoNTree(space, space, 0.1)
-##
-
-cfmm = CorrectionFactorMatrixMethod.PetrovGalerkinCFMM(op, space, space, tree)
-##
-x = rand(size(cfmm, 2))
-y = cfmm * x
-
-#A = assemble(op, space, space)
-using LinearAlgebra
-y2 = A * x
-norm(y - y2) / norm(y2)
diff --git a/test/h2trees.jl b/test/h2trees.jl
deleted file mode 100644
index bf428e4..0000000
--- a/test/h2trees.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-using H2Trees
-using CompScienceMeshes
-
-m = meshsphere(1.0, 0.1)
-tree = TwoNTree(vertices(m), 0.1; minvalues=60)
\ No newline at end of file
diff --git a/test/potentialmatrix.jl b/test/potentialmatrix.jl
deleted file mode 100644
index 0af350e..0000000
--- a/test/potentialmatrix.jl
+++ /dev/null
@@ -1,67 +0,0 @@
-using BEAST
-using CompScienceMeshes
-
-function sources(space::BEAST.Space, quadorder; dim=universedimension(space.geo))
-    elements, _, _ = assemblydata(space)
-    qp = BEAST.quadpoints(x -> refspace(space)(x), elements, (quadorder,))
-    points = zeros(Float64, length(qp) * length(qp[1, 1]), dim)
-
-    for (idx, pts) in enumerate(pts for element in qp for pts in element)
-        points[idx, 1:dim] = pts.point.cart[1:dim]
-    end
-
-    return points, qp
-end
-
-using OhMyThreads
-elements, _, _ = assemblydata(space)
-qp = BEAST.quadpoints(x -> refspace(space)(x), elements, (3,))
-OhMyThreads.@pmap for (idx, el) in enumerate(ts for el in qp for i in el)
-    println(idx)
-end
-
-ind = 1
-for el in qp
-    for pts in el
-        println(ind)
-        ind += 1
-    end
-end
-
-function potentialmatrix(qp::Q, X::BEAST.LagrangeBasis) where {Q}
-    rfspace = refspace(X)
-    _, tad, _ = assemblydata(X)
-    len = length(qp) * length(qp[1, 1]) * size(tad.data)[1] * size(tad.data)[2]
-
-    rc = ones(Int, len, 2)
-    vals = zeros(Float64, len)
-    sind = 1
-
-    for (ncell, cell) in enumerate(qp[1, :])
-        ind = (ncell - 1) * length(cell)
-        for (npoint, point) in enumerate(cell)
-            val = rfspace(point.point)
-            for localbasis in eachindex(val)
-                for data in tad.data[:, localbasis, ncell]
-                    if data[1] != 0 && ind + npoint != 0
-                        rc[sind, 1] = ind + npoint
-                        rc[sind, 2] = data[1]
-                        vals[sind] = val[localbasis].value * point.weight * data[2]
-                        sind += 1
-                    end
-                end
-            end
-        end
-    end
-
-    return rc, vals
-end
-
-##
-Γ = meshsphere(1.0, 0.08)
-
-op = Helmholtz3D.singlelayer()
-space = lagrangec0d1(Γ)
-typeof(space)
-
-points, qp = sources(space, 5)
diff --git a/test/runtests.jl b/test/runtests.jl
new file mode 100644
index 0000000..d2c5a50
--- /dev/null
+++ b/test/runtests.jl
@@ -0,0 +1,56 @@
+using Test
+using TestItems
+using TestItemRunner
+using CorrectionFactorMatrixMethod
+
+include("test_core.jl")
+
+using BEAST
+using CompScienceMeshes
+
+include("CFMMBEAST/test_CFMMBEAST.jl")
+
+using ExaFMMt
+
+include("CFMMExaFMMt/test_CFMMExaFMMt.jl")
+
+include("test_utils.jl")
+include("test_operators.jl")
+
+@testitem "Code quality (Aqua.jl)" begin
+    using Aqua
+    using CorrectionFactorMatrixMethod
+
+    Aqua.test_all(CorrectionFactorMatrixMethod)
+end
+
+@testitem "Code formatting (JuliaFormatter.jl)" begin
+    using JuliaFormatter
+    using CorrectionFactorMatrixMethod
+
+    @test JuliaFormatter.format(pkgdir(CorrectionFactorMatrixMethod); overwrite=false)
+end
+
+@testitem "Static analysis (JET.jl)" begin
+    using JET
+    using CorrectionFactorMatrixMethod
+
+    JET.test_package(
+        CorrectionFactorMatrixMethod; target_modules=(CorrectionFactorMatrixMethod,)
+    )
+end
+
+@testitem "Explicit imports (ExplicitImports.jl)" begin
+    using ExplicitImports
+    using CorrectionFactorMatrixMethod
+
+    @test ExplicitImports.check_no_stale_explicit_imports(CorrectionFactorMatrixMethod) ===
+        nothing
+    @test ExplicitImports.check_all_explicit_imports_via_owners(
+        CorrectionFactorMatrixMethod
+    ) === nothing
+    @test ExplicitImports.check_no_self_qualified_accesses(CorrectionFactorMatrixMethod) ===
+        nothing
+end
+
+@run_package_tests verbose = true
diff --git a/test/test_core.jl b/test/test_core.jl
new file mode 100644
index 0000000..d52d715
--- /dev/null
+++ b/test/test_core.jl
@@ -0,0 +1,111 @@
+using H2Trees
+using LinearMaps
+using StaticArrays
+
+const CFM = CorrectionFactorMatrixMethod
+
+struct TestCorrectedKernelMatrix{T} <: CFM.AbstractCorrectedKernelMatrix{T}
+    data::Matrix{T}
+end
+
+Base.size(matrix::TestCorrectedKernelMatrix, dim...) = size(matrix.data, dim...)
+
+function (matrix::TestCorrectedKernelMatrix)(block, testindices, trialindices)
+    block .= view(matrix.data, testindices, trialindices)
+    return nothing
+end
+
+struct TestAssembler
+    tfs::Vector{Int}
+    bfs::Vector{Int}
+end
+
+struct TestCorrectionOperator
+    farquadstrat
+end
+
+struct TestFMMFunctor <: CFM.FMMFunctor end
+
+CFM.scalartype(::TestCorrectionOperator) = Float64
+
+function CFM.AbstractCorrectedKernelMatrix(
+    operator::TestCorrectionOperator, testspace, trialspace; nearquadstrat, farquadstrat
+)
+    @assert farquadstrat === operator.farquadstrat
+    data = [10.0i + j for i in eachindex(testspace), j in eachindex(trialspace)]
+    return TestCorrectedKernelMatrix(data)
+end
+
+function (::TestFMMFunctor)(operator, testspace, trialspace; kwargs...)
+    return LinearMap(zeros(CFM.scalartype(operator), length(testspace), length(trialspace)))
+end
+
+@testset "Core" begin
+    quadrature = CFM.defaultfarquadstrat(nothing, nothing, nothing)
+    @test quadrature isa CFM.SafeDoubleNumQStrat
+    @test (quadrature.outer_rule, quadrature.inner_rule) == (3, 3)
+
+    data = [1.0 2.0 3.0; 4.0 5.0 6.0]
+    kernelmatrix = TestCorrectedKernelMatrix(data)
+    block = zeros(2, 2)
+    kernelmatrix(block, [2, 1], [1, 3])
+    @test eltype(kernelmatrix) == Float64
+    @test size(kernelmatrix) == size(data)
+    @test block == [4.0 6.0; 1.0 3.0]
+
+    assembler = TestAssembler(collect(1:3), collect(1:4))
+    beastmatrix = CFM.BEASTCorrectedKernelMatrix{ComplexF64}(assembler, assembler)
+    @test eltype(beastmatrix) == ComplexF64
+    @test size(beastmatrix) == (3, 4)
+
+    testspace = [SVector(x, 0.0, 0.0) for x in 0.0:3.0]
+    trialspace = [SVector(x, 0.0, 0.0) for x in 0.5:1.0:2.5]
+    tree = TwoNTree(testspace, trialspace, 0.25; testminvalues=1, trialminvalues=1)
+    farquadstrat = CFM.SafeDoubleNumQStrat(2, 3)
+    correctednears = CFM.assemblecorrectednears(
+        TestCorrectionOperator(farquadstrat),
+        testspace,
+        trialspace,
+        tree;
+        nearquadstrat=nothing,
+        farquadstrat=farquadstrat,
+        scheduler=CFM.SerialScheduler(),
+        isnear=(args...) -> true,
+    )
+    dense = [10.0i + j for i in eachindex(testspace), j in eachindex(trialspace)]
+    @test correctednears * ones(length(trialspace)) == dense * ones(length(trialspace))
+
+    assembled = CFMM.assemble(
+        TestCorrectionOperator(farquadstrat),
+        testspace,
+        trialspace;
+        fmmfunctor=TestFMMFunctor(),
+        minhalfsize=0.25,
+        testminvalues=1,
+        trialminvalues=1,
+        nearquadstrat=nothing,
+        farquadstrat=farquadstrat,
+        scheduler=CFM.SerialScheduler(),
+        isnear=(args...) -> true,
+    )
+    @test assembled * ones(length(trialspace)) == dense * ones(length(trialspace))
+    @test CFM.defaultminhalfsize(testspace, trialspace) == 1.5 / 2^10
+
+    options = CFM.FMMFunctor(; p=10, ncrit=64)
+    @test (options.p, options.ncrit) == (10, 64)
+    @test CFM.defaultminvalues(options) == 64
+    matrix = [2.0 1.0; -1.0 3.0]
+    vector = [1.0, 2.0]
+    @test CFM.FMM(matrix) * vector == matrix * vector
+    @test transpose(CFM.FMM(matrix, transpose(matrix))) * vector ==
+        transpose(matrix) * vector
+
+    far = LinearMap(matrix)
+    correction = LinearMap([0.5 -0.5; 1.0 0.0])
+    cfmm = CFM.PetrovGalerkinCFMM{Float64}(
+        correction, far, size(far), CFM.SerialScheduler()
+    )
+    @test cfmm * vector == (far + correction) * vector
+    @test transpose(cfmm) * vector == transpose(far + correction) * vector
+    @test adjoint(cfmm) * vector == adjoint(far + correction) * vector
+end
diff --git a/test/test_operators.jl b/test/test_operators.jl
new file mode 100644
index 0000000..b25cb95
--- /dev/null
+++ b/test/test_operators.jl
@@ -0,0 +1,25 @@
+@testset "Operators" begin
+    scalarproblems, vectorproblems = testproblems()
+    k = 2π / 10
+
+    for operator in (
+        Helmholtz3D.singlelayer(; wavenumber=k),
+        Helmholtz3D.doublelayer(; wavenumber=k),
+        Helmholtz3D.doublelayer_transposed(; wavenumber=k),
+        Helmholtz3D.hypersingular(; wavenumber=k),
+    )
+        testoperator(operator, scalarproblems)
+    end
+
+    for operator in
+        (Maxwell3D.singlelayer(; wavenumber=k), Maxwell3D.doublelayer(; wavenumber=k))
+        testoperator(operator, vectorproblems)
+    end
+
+    problem = last(vectorproblems)
+    operator = Maxwell3D.singlelayer(; wavenumber=k)
+    matrix = CFMM.assemble(operator, problem.testspace, problem.trialspace)
+    dense = assemble(operator, problem.testspace, problem.trialspace)
+    x = ones(ComplexF64, size(matrix, 2))
+    @test isapprox(matrix * x, dense * x; rtol=problem.tolerance)
+end
diff --git a/test/test_utils.jl b/test/test_utils.jl
new file mode 100644
index 0000000..940347f
--- /dev/null
+++ b/test/test_utils.jl
@@ -0,0 +1,81 @@
+using H2Trees
+using LinearAlgebra
+using StaticArrays
+
+function testoperator(operator, problems)
+    for (; testspace, trialspace, tree, tolerance) in problems
+        dense = assemble(operator, testspace, trialspace)
+        cfmm = CFMM.assemble(operator, testspace, trialspace; tree)
+
+        for matrix in (identity, transpose, adjoint), T in (ComplexF64, ComplexF32)
+            n = size(matrix(cfmm), 2)
+            x = T.(complex.((1:n) ./ n, (n:-1:1) ./ n))
+            y = matrix(cfmm) * x
+
+            @test eltype(y) == promote_type(eltype(cfmm), T)
+            @test isapprox(y, matrix(dense) * x; rtol=tolerance)
+        end
+    end
+end
+
+function testproblems()
+    source = meshrectangle(1.0, 1.0, 0.3)
+    scalarspace = lagrangec0d1(source)
+    vectorspace = raviartthomas(source)
+    scalarproblems = NamedTuple[]
+    vectorproblems = NamedTuple[]
+
+    for (offset, tolerance) in
+        ((SVector(3.0, 0.0, 1.0), 2.0e-4), (SVector(0.2, 0.0, 0.2), 2.0e-3))
+        target = translate(source, offset)
+        testscalar = lagrangec0d1(target)
+        testvector = raviartthomas(target)
+        push!(
+            scalarproblems,
+            (
+                testspace=testscalar,
+                trialspace=scalarspace,
+                tree=TwoNTree(
+                    testscalar, scalarspace, 0.1; testminvalues=50, trialminvalues=50
+                ),
+                tolerance=tolerance,
+            ),
+        )
+        push!(
+            vectorproblems,
+            (
+                testspace=testvector,
+                trialspace=vectorspace,
+                tree=TwoNTree(
+                    testvector, vectorspace, 0.1; testminvalues=50, trialminvalues=50
+                ),
+                tolerance=tolerance,
+            ),
+        )
+    end
+
+    push!(
+        scalarproblems,
+        (
+            testspace=scalarspace,
+            trialspace=scalarspace,
+            tree=TwoNTree(
+                scalarspace, scalarspace, 0.1; testminvalues=50, trialminvalues=50
+            ),
+            tolerance=2.0e-4,
+        ),
+    )
+    push!(
+        vectorproblems,
+        (
+            testspace=vectorspace,
+            trialspace=vectorspace,
+            tree=TwoNTree(
+                vectorspace, vectorspace, 0.1; testminvalues=50, trialminvalues=50
+            ),
+            tolerance=2.0e-4,
+        ),
+    )
+
+    return scalarproblems, vectorproblems
+end