diff --git a/docs/src/examples/helmholtz_scattering.jl b/docs/src/examples/helmholtz_scattering.jl
index b262b3ac..25ee13b6 100644
--- a/docs/src/examples/helmholtz_scattering.jl
+++ b/docs/src/examples/helmholtz_scattering.jl
@@ -56,7 +56,7 @@ Pkg.activate(docsdir)                 #src
 # and the *Sommerfeld radiation condition* at infinity
 
 # ```math
-#     \lim_{|\boldsymbol{x}| \to \infty} \|\boldsymbol{x}|^{(d-1)/2} \left( \frac{\partial u}{\partial |\boldsymbol{x}|} - i k u \right) = 0.
+#     \lim_{|\boldsymbol{x}| \to \infty} |\boldsymbol{x}|^{(d-1)/2} \left( \frac{\partial u}{\partial |\boldsymbol{x}|} - i k u \right) = 0.
 # ```
 
 # Here ``g`` is a (given) boundary datum, and ``k`` is the constant wavenumber.
diff --git a/docs/src/index.md b/docs/src/index.md
index 299d1d01..6b391b4a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,15 +37,20 @@ the latest stable version of Julia, although `Inti.jl` should work with
 
 ## Installing Inti.jl
 
-Inti.jl is not yet registered in the Julia General registry. You can install it
-by launching a Julia REPL and typing the following command:
+Inti.jl is registered in the Julia General registry and can be installed by 
+launching a Julia REPL and typing the following command:
+
+```julia
+]add Inti
+```
+
+Alternatively, you can install the latest version of Inti.jl from the `main` branch using:
 
 ```julia
 using Pkg; Pkg.add(;url = "https://github.com/IntegralEquations/Inti.jl", rev = "main")
 ```
 
-This will download and install the latest version of Inti.jl from the `main`
-branch. Change `rev` if you need a different branch or a specific commit hash.
+Change `rev` if you need a different branch or a specific commit hash.
 
 ## Installing weak dependencies
 
diff --git a/docs/src/tutorials/compression_methods.md b/docs/src/tutorials/compression_methods.md
index 94a79d7c..e2daedaf 100644
--- a/docs/src/tutorials/compression_methods.md
+++ b/docs/src/tutorials/compression_methods.md
@@ -23,7 +23,7 @@ of the operator. The following methods are available:
   [`HMatrices`](https://github.com/IntegralEquations/HMatrices.jl) library.
 - [`assemble_fmm`](@ref): return a `LinearMap` object that represents the
   operator using the fast multipole method. This method is powered by the
-  [`FMM2D`](https://github.com/flatironinstitute/fmm2d/) and
+  [`FMM2D`](https://github.com/flatironinstitute/fmm2d/), [`FMMLIB2D`](https://github.com/zgimbutas/fmmlib2d) and
   [`FMM3D`](https://fmm3d.readthedocs.io) libraries, and is only available for
   certain kernels.
 
@@ -98,7 +98,7 @@ println("Forward map error: $er")
 
 Note that `HMatrices` are said to be *kernel-independent*, meaning that they
 efficiently compress a wide range of integral operators provided they satisfy a
-certain asymptotically smooth criterion (see e.g. [bebendorf2008hierarchical,
+certain asymptotic smoothness criterion (see e.g. [bebendorf2008hierarchical,
 hackbusch2015hierarchical](@cite)).
 
 The `HMatrix` object can be used to solve linear systems, both iteratively
diff --git a/docs/src/tutorials/geo_and_meshes.md b/docs/src/tutorials/geo_and_meshes.md
index a8999fba..49b7e2a0 100644
--- a/docs/src/tutorials/geo_and_meshes.md
+++ b/docs/src/tutorials/geo_and_meshes.md
@@ -28,7 +28,7 @@ function that maps points from a [`ReferenceShape`](@ref) to the physical space.
 In most applications involving complex three-dimensional surfaces, an external
 meshing software is used to generate a mesh, and the mesh is imported using the
 `import_mesh` function (which relies on [Gmsh](https://gmsh.info)). The entities
-can then the extracted from the mesh based on e.g. their dimension or label.
+can then be extracted from the mesh based on e.g. their dimension or label.
 Here is an example of how to import a mesh from a file:
 
 ```@example geo-and-meshes
@@ -54,7 +54,7 @@ nothing # hide
 ```
 
 You can filter entities satisfying a certain condition, e.g., entities of a
-given dimension of containing a certain label, in order to construct a domain:
+given dimension or containing a certain label, in order to construct a domain:
 
 ```@example geo-and-meshes
 filter = e -> Inti.geometric_dimension(e) == 3
@@ -213,7 +213,7 @@ including the boundary segments:
 Inti.entities(msh)
 ```
 
-This allows you into the `msh` object to extract e.g. the boundary mesh:
+This allows you to probe the `msh` object to extract e.g. the boundary mesh:
 
 ```@example geo-and-meshes
 viz(msh[Inti.boundary(Ω)]; color = :red)
diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md
index c612d22d..8e0570d0 100644
--- a/docs/src/tutorials/getting_started.md
+++ b/docs/src/tutorials/getting_started.md
@@ -101,7 +101,7 @@ Q[1]
 
 In the constructor above we specified a quadrature order of 5, and Inti.jl
 internally picked a [`ReferenceQuadrature`](@ref) suitable for the specified
-order; for a finer control, you can also specify a quadrature rule directly.
+order; for finer control, you can also specify a quadrature rule directly.
 
 ## Integral operators
 
@@ -167,9 +167,10 @@ equation. For that, we need to provide the boundary data ``g``.
 We are interested in the scattered field ``u`` produced by an incident plane
 wave ``u_i = e^{i k \boldsymbol{d} \cdot \boldsymbol{x}}``, where
 ``\boldsymbol{d}`` is a unit vector denoting the direction of the plane wave.
-Assuming that the total field ``u_t = u_i + u`` satisfies a homogenous Neumann
-condition on ``\Gamma``, and that the scattered field ``u`` satisfies the
-Sommerfeld radiation condition, we can write the boundary condition as:
+Assuming that the total field ``u_t = u_i + u``, unique up to a constant,
+satisfies a homogenous Neumann condition on ``\Gamma``, and that the scattered
+field ``u`` satisfies the Sommerfeld radiation condition, we can write the
+boundary condition as:
 
 ```math
     \partial_\nu u = -\partial_\nu u_i, \quad \boldsymbol{x} \in \Gamma.
diff --git a/docs/src/tutorials/layer_potentials.md b/docs/src/tutorials/layer_potentials.md
index cce58765..dcf1d302 100644
--- a/docs/src/tutorials/layer_potentials.md
+++ b/docs/src/tutorials/layer_potentials.md
@@ -26,7 +26,7 @@ where ``K`` is the kernel of the operator, ``\Gamma`` is the source's boundary,
 ``\boldsymbol{r} \not \in \Gamma`` is a target point, and ``\sigma`` is the
 source density.
 
-Here is a simple example of how to create a kernel representing Laplace's
+Here is a simple example of how to create a kernel representing a Laplace
 double-layer potential:
 
 ```@example layer_potentials
@@ -44,8 +44,8 @@ Q = Inti.Quadrature(Γ; meshsize = 0.1, qorder = 5)
 𝒮 = Inti.IntegralPotential(K, Q)
 ```
 
-If you have a source density ``\sigma``, defined on the quadrature nodes of
-``\Gamma``, you can create a function that evaluates the layer potential at an
+If we have a source density ``\sigma``, defined on the quadrature nodes of
+``\Gamma``, we can create a function that evaluates the layer potential at an
 arbitrary point:
 
 ```@example layer_potentials
@@ -141,9 +141,9 @@ representation holds:
 u(\boldsymbol{r}) = \mathcal{S}[\gamma_1 u](\boldsymbol{r}) - \mathcal{D}[\gamma_0 u](\boldsymbol{r}), \quad \boldsymbol{r} \in \Omega
 ```
 
-where ``\gamma_0 u`` and ``\gamma_1 u`` are the Dirichlet and Neumann traces of
-``u``, and ``\mathcal{S}`` and ``\mathcal{D}`` are the single and double layer
-potentials over ``\Gamma := \partial \Omega``.
+where ``\gamma_0 u`` and ``\gamma_1 u`` are the respective Dirichlet and
+Neumann traces of ``u``, and ``\mathcal{S}`` and ``\mathcal{D}`` are the respective
+single and double layer potentials over ``\Gamma := \partial \Omega``.
 
 Let's compare next the exact solution with the layer potential evaluation, based
 on a quadrature of ``\Gamma``:
@@ -185,10 +185,11 @@ recommended:
 2. When you wish to evaluate the layer potential at many target points and take
    advantage of an acceleration routine.
 
-In such cases, it is recommended to use the `single_double_layer` function, (or
-to directly assemble an `IntegralOperator`) with a correction method. Here is an
-example of how to use the FMM acceleration with a near-field correction to
-evaluate the layer potentials::
+In such cases, it is recommended to use the `single_double_layer` function
+(alternately, one can directly assemble an `IntegralOperator`) with a
+correction and/or compression method as appropriate. Here is an example of how
+to use the FMM acceleration with a near-field correction to evaluate the layer
+potentials::
 
 ```@example layer_potentials
 using FMM2D
diff --git a/test/convergence/dim_convergence.jl b/test/convergence/dim_convergence.jl
index a9c13cb3..25adf384 100644
--- a/test/convergence/dim_convergence.jl
+++ b/test/convergence/dim_convergence.jl
@@ -11,6 +11,7 @@ atol = 0
 rtol = 1e-8
 t = :exterior
 σ = t == :interior ? 1 / 2 : -1 / 2
+# For N = 3 one needs to use compression, e.g. `compression = :fmm`, for the operators below
 N = 2
 pde = Inti.Laplace(; dim = N)
 # pde = Inti.Helmholtz(; dim = N, k = 2π)
@@ -35,12 +36,16 @@ for h in hh
     gmsh.option.setNumber("General.Verbosity", 2)
     gmsh.model.mesh.setOrder(2)
     Inti.clear_entities!()
-    gmsh.model.occ.addDisk(center[1], center[2], 0, 2 * radius, radius)
+    if N == 2
+        gmsh.model.occ.addDisk(center[1], center[2], 0, 2 * radius, radius)
+    else
+        gmsh.model.occ.addSphere(center[1], center[2], center[3], radius)
+    end
     gmsh.model.occ.synchronize()
     gmsh.model.mesh.generate(2)
-    M = Inti.import_mesh(; dim = 2)
+    M = Inti.import_mesh(; dim = N)
     gmsh.finalize()
-    Γ = Inti.Domain(e -> Inti.geometric_dimension(e) == 1, Inti.entities(M))
+    Γ = Inti.Domain(e -> Inti.geometric_dimension(e) == N - 1, Inti.entities(M))
     Q = Inti.Quadrature(view(M, Γ); qorder)
     @show Q
     xs = if t == :interior
@@ -83,10 +88,16 @@ end
 order = n + 1
 fig = Figure()
 ax = Axis(fig[1, 1]; xscale = log10, yscale = log10, xlabel = "h", ylabel = "error")
-scatterlines!(ax, hh, ee0; m = :x, label = "nocorrection")
-scatterlines!(ax, hh, ee1; m = :x, label = "dim correction")
+scatterlines!(ax, hh, ee0; marker = :x, label = "nocorrection")
+scatterlines!(ax, hh, ee1; marker = :x, label = "dim correction")
 ref = hh .^ order
 iref = length(ref)
-lines!(ax, hh, ee1[iref] / ref[iref] * ref; label = L"\mathcal{O}(h^%$order)", ls = :dash)
+lines!(
+    ax,
+    hh,
+    ee1[iref] / ref[iref] * ref;
+    label = L"\mathcal{O}(h^%$order)",
+    linestyle = :dash,
+)
 axislegend(ax)
 fig