JoshuaTetzner · JoshuaTetzner · Jul 2, 2026 · Jul 1, 2026 · Jul 1, 2026 · Jul 1, 2026
diff --git a/.gitignore b/.gitignore
@@ -10,3 +10,4 @@ lcov.info
 /docs/Manifest.toml
 /docs/build/
 /docs/site/
+/examples/*.html
diff --git a/README.md b/README.md
@@ -22,21 +22,18 @@ extensions for
 [ExaFMMt.jl](https://github.com/JoshuaTetzner/ExaFMMt.jl), which are loaded
 automatically once these packages are available.
 
+Planned extensions include interfaces to additional FMM backends — a
+native Julia implementation and GPU-accelerated variants — as well as support
+for two-dimensional problems.
+
 ## 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 FMM [1]:
 
-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/).
+$$\boldsymbol{A}\boldsymbol{x} \approx \boldsymbol{A}_\mathrm{FMM}\,\boldsymbol{x} + \left(\boldsymbol{A}_\mathrm{near} - \boldsymbol{A}_\mathrm{FMM,near}\right)\boldsymbol{x}.$$
 
 ## Installation
 
@@ -63,18 +60,10 @@ mesh = meshsphere(1.0, 0.4)
 space = raviartthomas(mesh)
 operator = Maxwell3D.singlelayer(; wavenumber=1.0)
 
-matrix = CFMM.assemble(operator, space)
+matrix = CFMM.assemble(operator, space, 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/make.jl b/docs/make.jl
@@ -24,6 +24,7 @@ makedocs(;
             ["General Usage" => "manual/manual.md", "Examples" => "manual/examples.md"],
         "Details" => [
             "Correction-Factor Method" => "details/method.md",
+            "Supported Operators" => "details/operators.md",
             "Package Extensions" => "details/extensions.md",
         ],
         "Contributing" => "contributing.md",

diff --git a/docs/render_examples.jl b/docs/render_examples.jl
@@ -0,0 +1,28 @@
+# Regenerate the pre-rendered plot assets embedded in the documentation.
+# Run from the package root:
+#
+#   julia --startup-file=no docs/render_examples.jl
+#
+# Commit the generated files in docs/src/assets/examples/ afterwards.
+
+using Pkg: Pkg
+Pkg.activate(joinpath(@__DIR__, "..", "examples"))
+Pkg.develop(; path=joinpath(@__DIR__, ".."))
+Pkg.instantiate()
+
+using BEAST
+using CompScienceMeshes
+using CorrectionFactorMatrixMethod
+using ExaFMMt
+using Krylov
+using LinearAlgebra
+using PlotlyJS
+
+outdir = joinpath(@__DIR__, "src", "assets", "examples")
+mkpath(outdir)
+ENV["CFMM_OUTPUT_DIR"] = outdir
+
+include(joinpath(@__DIR__, "..", "examples", "efie.jl"))
+include(joinpath(@__DIR__, "..", "examples", "mfie.jl"))
+
+@info "Plots written to $outdir"
diff --git a/docs/src/assets/examples/efie_results.html b/docs/src/assets/examples/efie_results.html
diff --git a/docs/src/assets/examples/mfie_results.html b/docs/src/assets/examples/mfie_results.html
diff --git a/docs/src/contributing.md b/docs/src/contributing.md
@@ -1,7 +1,95 @@
 # 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.
+## What is currently implemented
+
+CorrectionFactorMatrixMethod.jl provides a matrix-free BEM framework for
+**three-dimensional** problems. The following is in place:
+
+- **Core machinery**: correction-factor assembly driven by an `H2Trees.BlockTree`,
+  sparse near-field correction via `BlockSparseMatrices.jl`, and a
+  `LinearMaps.jl`-compatible `PetrovGalerkinCFMM` operator supporting forward,
+  transpose, and adjoint products.
+- **BEAST extension** (`CFMMBEAST`): quadrature, collocation, and corrected-kernel
+  assembly for all BEAST `IntegralOperator` types; specialisations for
+  Helmholtz 3D (V, K, K', W) with Lagrange C0D1 basis functions and Maxwell 3D
+  (T/EFIE, K/MFIE) with Raviart–Thomas basis functions.
+- **ExaFMMt extension** (`CFMMExaFMMt`): operator-parameter mapping to
+  `LaplaceFMMOptions`, `ModifiedHelmholtzFMMOptions`, and `HelmholtzFMMOptions`;
+  Helmholtz operators with zero wavenumber automatically reduce to the Laplace
+  path.
+
+For a complete list of tested operator–basis combinations see the
+[Supported Operators](@ref) page.
+
+## What could be added
+
+Contributions are welcome in any of the following areas:
+
+### Additional FMM backends
+
+The `FMMFunctor` interface is designed to be extended. Adding a new FMM
+library requires implementing:
+
+- `CorrectionFactorMatrixMethod.setup(spoints, tpoints, options)` — set up the
+  FMM from source and target point matrices.
+- `CorrectionFactorMatrixMethod.fmmresult(fmm, x)` — evaluate the FMM
+  map for a given input vector.
+
+Potential targets include:
+
+- **A native Julia FMM**: A pure-Julia implementation would remove the binary
+  dependency on ExaFMMt and simplify installation, particularly on HPC
+  systems.
+- **FMMLIB2D / FMMLIB3D**: The Fortran-based FMM libraries wrapped for Julia
+  could serve as an alternative backend.
+- **GPU-accelerated FMMs**: An extension that drives a GPU FMM (e.g. via CUDA)
+  would be a natural fit for the GPU near-field assembler already used for
+  BEAST.
+
+### Two-dimensional operator support
+
+The current collocation and operator implementations target 3D surfaces only.
+Extending to 2D would require:
+
+- A 2D `sources` / `potentials` implementation for 1D boundary curves.
+- Operators for the 2D Helmholtz and Laplace kernels
+  (``G(\boldsymbol{x},\boldsymbol{y}) = \tfrac{i}{4} H_0^{(1)}(k|\boldsymbol{x}-\boldsymbol{y}|)``).
+- A 2D-capable FMM backend (e.g. FMMLIB2D).
+
+### Additional BEM library integrations
+
+The `CFMMBEAST` extension can serve as a blueprint for integrating other
+Julia BEM libraries. Any library that can provide quadrature data and
+block assemblers in the same style can be added as a new extension.
+
+### Further operator–basis combinations
+
+Some operator–basis pairings (e.g. Helmholtz operators with
+Raviart–Thomas basis functions, or Maxwell operators with Lagrange
+basis functions for mixed formulations) are architecturally supported but not
+yet tested. Contributions that add and verify such combinations are welcome.
+
+## Regenerating documentation plots
+
+The interactive plots on the Examples page are pre-rendered HTML files stored
+in `docs/src/assets/examples/` and committed to the repository. The CI build
+only serves these static files; it does not re-run the examples.
+
+After changing any file in `examples/`, regenerate the plots from the
+repository root:
+
+```sh
+julia --startup-file=no docs/render_examples.jl
+```
+
+Then commit the updated HTML files alongside your example changes. Running the
+test suite locally will emit a `@warn` if the HTML files are older than the
+example scripts, so you will not accidentally forget this step.
+
+## Development workflow
+
+Development targets the `dev` branch. Open a pull request from `dev` to
+`main` for stable releases and merge with a merge commit.
 
 Run the package checks locally with:
 
@@ -20,4 +108,5 @@ include("docs/make.jl")
 ```
 
 Code must pass the formatter, Aqua, JET, and ExplicitImports checks included in
-the test suite.
+the test suite. The formatter is configured with the `blue` style at a column
+margin of 92.
diff --git a/docs/src/details/method.md b/docs/src/details/method.md
@@ -6,8 +6,8 @@ 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.
+\boldsymbol{A}\boldsymbol{x} \approx \boldsymbol{A}_{\mathrm{FMM}}\boldsymbol{x} +
+\left(\boldsymbol{A}_{\mathrm{near}} - \boldsymbol{A}_{\mathrm{FMM,near}}\right)\boldsymbol{x}.
 ```
 
 The sparse correction blocks are selected from an `H2Trees.BlockTree`.

diff --git a/docs/src/details/operators.md b/docs/src/details/operators.md
@@ -0,0 +1,69 @@
+# Supported Operators
+
+All operators currently implemented target **three-dimensional** problems.
+Two-dimensional support is a planned future extension.
+
+The CFMM wrapper is constructed via [`CFMM.assemble`](@ref) (or the lower-level
+[`PetrovGalerkinCFMM`](@ref)) once the `CFMMBEAST` and `CFMMExaFMMt` extensions
+are loaded.
+
+## Helmholtz 3D (scalar)
+
+These operators act on scalar basis functions.
+The current implementation supports **Lagrange C0D1** (`lagrangec0d1`) basis
+functions on triangulated surfaces.
+
+The Helmholtz Green's function is
+
+```math
+G(\boldsymbol{x},\boldsymbol{y}) = \frac{e^{ik|\boldsymbol{x}-\boldsymbol{y}|}}{4\pi|\boldsymbol{x}-\boldsymbol{y}|},
+\quad k \in \mathbb{C}.
+```
+
+Setting `wavenumber = 0` recovers the Laplace Green's function
+``G_0(\boldsymbol{x},\boldsymbol{y}) = 1/(4\pi|\boldsymbol{x}-\boldsymbol{y}|)``; the ExaFMMt backend
+automatically switches to `LaplaceFMMOptions` in that case.
+
+| Operator | BEAST constructor | Mathematical form |
+|:---------|:------------------|:------------------|
+| Single-layer (V) | `Helmholtz3D.singlelayer(; wavenumber)` | ``(V u)(\boldsymbol{x}) = \int_\Gamma G(\boldsymbol{x},\boldsymbol{y})\,u(\boldsymbol{y})\,\mathrm{d}S(\boldsymbol{y})`` |
+| Double-layer (K) | `Helmholtz3D.doublelayer(; wavenumber)` | ``(K u)(\boldsymbol{x}) = \int_\Gamma \partial_{\boldsymbol{n}(\boldsymbol{y})} G(\boldsymbol{x},\boldsymbol{y})\,u(\boldsymbol{y})\,\mathrm{d}S(\boldsymbol{y})`` |
+| Transposed double-layer (K') | `Helmholtz3D.doublelayer_transposed(; wavenumber)` | ``(K' u)(\boldsymbol{x}) = \int_\Gamma \partial_{\boldsymbol{n}(\boldsymbol{x})} G(\boldsymbol{x},\boldsymbol{y})\,u(\boldsymbol{y})\,\mathrm{d}S(\boldsymbol{y})`` |
+| Hypersingular (W) | `Helmholtz3D.hypersingular(; wavenumber)` | ``(W u)(\boldsymbol{x}) = -\partial_{\boldsymbol{n}(\boldsymbol{x})} \int_\Gamma \partial_{\boldsymbol{n}(\boldsymbol{y})} G(\boldsymbol{x},\boldsymbol{y})\,u(\boldsymbol{y})\,\mathrm{d}S(\boldsymbol{y})`` |
+
+All four operators are tested with `lagrangec0d1` basis functions (Petrov–Galerkin,
+square systems). Other scalar basis function types may work but are not
+currently verified.
+
+## Maxwell 3D (vector)
+
+These operators act on vector basis functions.
+The current implementation supports **Raviart–Thomas** (`raviartthomas`) basis
+functions on triangulated surfaces.
+
+| Operator | BEAST constructor | Mathematical form |
+|:---------|:------------------|:------------------|
+| Single-layer / EFIE (T) | `Maxwell3D.singlelayer(; wavenumber)` | ``(\mathcal{T}\boldsymbol{u})(\boldsymbol{x}) = \alpha\!\int_\Gamma G\,\boldsymbol{u}\,\mathrm{d}S + \beta\,\nabla\!\int_\Gamma G\,\nabla_\Gamma\!\cdot\boldsymbol{u}\,\mathrm{d}S`` |
+| Double-layer / MFIE (K) | `Maxwell3D.doublelayer(; wavenumber)` | ``(\mathcal{K}\boldsymbol{u})(\boldsymbol{x}) = \int_\Gamma \nabla G(\boldsymbol{x},\boldsymbol{y})\times\boldsymbol{u}(\boldsymbol{y})\,\mathrm{d}S(\boldsymbol{y})`` |
+
+Both operators are tested with `raviartthomas` basis functions (Petrov–Galerkin,
+square systems).
+
+## Tested basis-function / operator combinations
+
+The table below summarises which combinations are covered by the automated
+test suite. A ✓ means the forward product, the transpose, and the adjoint
+are all verified against a dense BEAST assembly.
+
+| Operator | `lagrangec0d1` | `raviartthomas` |
+|:---------|:--------------:|:---------------:|
+| `Helmholtz3D.singlelayer` | ✓ | — |
+| `Helmholtz3D.doublelayer` | ✓ | — |
+| `Helmholtz3D.doublelayer_transposed` | ✓ | — |
+| `Helmholtz3D.hypersingular` | ✓ | — |
+| `Maxwell3D.singlelayer` | — | ✓ |
+| `Maxwell3D.doublelayer` | — | ✓ |
+
+Combinations marked — are not supported by the current operator
+implementation (the scalar Helmholtz operators require scalar basis functions;
+the vector Maxwell operators require vector basis functions).
diff --git a/docs/src/manual/examples.md b/docs/src/manual/examples.md
@@ -36,16 +36,26 @@ space = raviartthomas(mesh)
 wavenumber = 1.0
 operator = Maxwell3D.singlelayer(; wavenumber)
 excitation = Maxwell3D.planewave(;
-    direction=ẑ, polarization=x̂, wavenumber
+    direction=ẑ, polarization=x̂, wavenumber
 )
 rhs = assemble((n × excitation) × n, space)
 
-matrix = CFMM.assemble(operator, space)
+matrix = CFMM.assemble(operator, space, space)
 current, stats = Krylov.gmres(
     matrix, rhs; rtol=1.0e-4, itmax=200, history=true, verbose=1
 )
 ```
 
+Far-field pattern (top left), scattered electric field magnitude in the ``yz``
+plane (top right), and surface-current magnitude (bottom):
+
+```@raw html
+<iframe src="../../assets/examples/efie_results.html"
+        style="width:100%;height:600px;border:none;"
+        loading="lazy">
+</iframe>
+```
+
 The complete runnable example, including the plots, is in
 `examples/efie.jl`.
 
@@ -67,7 +77,7 @@ wavenumber = frequency * √(permittivity * permeability)
 
 operator = Maxwell3D.doublelayer(; wavenumber)
 excitation = Maxwell3D.planewave(;
-    direction=ẑ, polarization=x̂, wavenumber
+    direction=ẑ, polarization=x̂, wavenumber
 )
 magneticfield = -1 / (im * permeability * frequency) * curl(excitation)
 rhs = assemble((n × magneticfield) × n, testspace)
@@ -80,5 +90,15 @@ current, stats = Krylov.gmres(
 )
 ```
 
+Far-field pattern (top left), scattered electric field magnitude in the ``yz``
+plane (top right), and surface-current magnitude (bottom):
+
+```@raw html
+<iframe src="../../assets/examples/mfie_results.html"
+        style="width:100%;height:600px;border:none;"
+        loading="lazy">
+</iframe>
+```
+
 The complete runnable example, including the plots, is in
 `examples/mfie.jl`.