Improve documentation (#165)

* Remove current docs index.md

* New index.md as symlink to README

* Update example/demo

* Rearrange and relabel the Support Matrix, update notes

* Change naming of Support Status

* Overhaul usage into new Tutorial page

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <[email protected]>

* Add a link to juliadocs for do syntax

* Convert to example blocks

Co-authored-by: Joshua Lampert <[email protected]>

* Add Meshes and Unitful to docs Project

* Convert another block

* Try a jldoctest

* Add repl prompts

* Revert to static code block

* Reorganize files

* Add a placeholder page for integration rules

* Add some draft contents

* Add some more notes

* Remove unneeded comment

* Add a Tips page

* Add links to sections

* fix typo

* Apply suggestion

Co-authored-by: Joshua Lampert <[email protected]>

* Change page order

* Clarify that atol tip affects both adaptive solvers

* Add more notes to tip

* Update language

* Add a tip for AutoEnzyme

* Update docs/src/tips.md

Co-authored-by: Joshua Lampert <[email protected]>

* Convert code blocks to repl blocks

* Add compats

* Update docs/Project.toml

Co-authored-by: Joshua Lampert <[email protected]>

---------

Co-authored-by: Joshua Lampert <[email protected]>

Author: mikeingold

README.md

@@ -30,7 +30,7 @@ Performs a numerical integration of some integrand function `f` over the domain
 ```julia
 integral(f, geometry, rule)
 ```
-Performs a numerical integration of some integrand function `f` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. The integrand function can be anything callable with a method `f(::Meshes.Point)`.  
+Performs a numerical integration of some integrand function `f` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. The integrand function can be anything callable with a method `f(::Meshes.Point)`.
 
 Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics.
 
@@ -46,24 +46,30 @@ Alias functions are provided for convenience. These are simply wrappers for `int
 - `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
 - `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)
 
-# Demo
+### Example
 
 ```julia
 using Meshes
 using MeshIntegrals
 using Unitful
 
-# Define a path that approximates a sine-wave on the xy-plane
-mypath = BezierCurve(
-    [Point(t*u"m", sin(t)*u"m", 0.0u"m") for t in range(-pi, pi, length=361)]
-)
+# Define a Bezier curve whose path approximates a sine-wave on the xy-plane
+N = 361
+curve = BezierCurve([Point(t*u"m", sin(t)*u"m", 0.0u"m") for t in range(-π, π, length=N)])
 
-# Map f(::Point) -> f(x, y, z) in unitless coordinates
-f(p::Meshes.Point) = f(ustrip(to(p))...)
+# Integrand function that outputs in units of Ohms/meter
+function f(p::Point)
+    x, y, z = to(p)
+    return (1 / sqrt(1 + cos(x / u"m")^2)) * u"Ω/m"
+end
 
-# Integrand function in units of Ohms/meter
-f(x, y, z) = (1 / sqrt(1 + cos(x)^2)) * u"Ω/m"
+# Use recommended defaults
+integral(f, curve) # -> Approximately 2π Ω
 
-integral(f, mypath)
-# -> Approximately 2*Pi Ω
+# Using aliases
+lineintegral(f, curve) # -> Approximately 2π Ω
+surfaceintegral(f, curve) # -> throws ArgumentError
+
+# Specifying an integration rule and settings (loosened absolute tolerance)
+integral(f, curve, GaussKronrod(atol = 1e-4u"Ω")) # -> Approximately (2π ± 1e-4) Ω
 ```

docs/Project.toml

@@ -1,3 +1,13 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 MeshIntegrals = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+
+[compat]
+BenchmarkTools = "1"
+Documenter = "1"
+Meshes = "0.51.20, 0.52"
+Unitful = "1.19"
+julia = "1.9"

docs/make.jl

@@ -11,13 +11,15 @@ makedocs(
     pages = [
         "Home" => [
             "About" => "index.md",
-            "Support Matrix" => "supportmatrix.md",
-            "Example Usage" => "usage.md"
+            "Tutorial" => "tutorial.md",
+            "Integration Rules" => "integration_rules.md",
+            "Support Status" => "support.md",
+            "Tips" => "tips.md"
         ],
         "Developer Notes" => [
-            "Changelog" => "CHANGELOG.md",
-            "How it Works" => "how_it_works.md",
-            "Specializations" => "specializations.md"
+            "Changelog" => "developer/CHANGELOG.md",
+            "How it Works" => "developer/how_it_works.md",
+            "Specializations" => "developer/specializations.md"
         ],
         "Public API" => "api.md"
     ]

docs/src/CHANGELOG.md

@@ -0,0 +1 @@
+../../../CHANGELOG.md

docs/src/developer/CHANGELOG.md

@@ -1,4 +1,4 @@
-# How it Works (By Example)
+# [How it Works](@id howitworks)
 
 ## Example Problem
 

docs/src/how_it_works.md docs/src/developer/how_it_works.mddocs/src/how_it_works.md renamed to docs/src/developer/how_it_works.md

@@ -0,0 +1 @@
+../../README.md

docs/src/specializations.md docs/src/developer/specializations.mddocs/src/specializations.md renamed to docs/src/developer/specializations.md

@@ -0,0 +1,55 @@
+# Integration Rules
+
+An **integration rule** is a method or algorithm used to numerically calculate the value of an integral. When an integral is calculated using the two-argument form `integral(f, geometry)`, **MeshIntegrals.jl** will automatically select an integration rule. Default rules are generally well-behaved, but may not be optimal for every problem.
+
+**MeshIntegrals.jl** defines an abstract type `IntegrationRule` with sub-types representing the various integration rules supported by this package. These rules can be specified when calculating an integral using the three-argument form `integral(f, geometry, rule)`. Currently, the following rule types are implemented:
+- [`GaussKronrod`](@ref gausskronrod) for adaptive Gauss-Kronrod quadrature rules
+- [`GaussLegendre`](@ref gausslegendre) for Gauss-Legendre quadrature rules
+- [`HAdaptiveCubature`](@ref hadaptivecubature) for h-adaptive cubature rules
+
+## [Gauss-Kronrod](@id gausskronrod)
+
+The `GaussKronrod` type is used to specify an adaptive Gauss-Kronrod quadrature rule, as implemented by [QuadGK.jl](https://github.com/JuliaMath/QuadGK.jl).
+
+```julia
+rule = GaussKronrod()
+```
+
+All standard `QuadGK.quadgk` keyword-argument options are supported. These can be specified when constructing the rule, where the `kwargs` in `GaussKronrod(kwargs...)` is equivalent to `quadgk(f, a, b; kwargs...)`, e.g.:
+```julia
+rule = GaussKronrod(order = 5, rtol = 1e-4)
+```
+
+## [Gauss-Legendre](@id gausslegendre)
+
+The `GaussLegendre` type is used to specify a [Gauss-Legendre quadrature](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature) rule. Gauss-Legendre quadrature rules of order $N$ are used to approximate definite integrals by sampling the integrand on a fixed grid with corresponding nodes $x_i$ and weights $w_i$.
+```math
+\int_{-1}^1 f(x) ~\text{d}x \approx \sum_{i=1}^N w_i \, f(x_i)
+```
+
+These nodes and weights are purely a function of $N$ as they are derived from the roots of $N$-th order Legendre polynomials. This has the effect of providing exact solutions for integrand functions $f$ that can be represented as degree $2N-1$ polynomials.
+
+This integration process can also be extended into multiple dimensions, for example:
+```math
+\int_{-1}^1 \int_{-1}^1 \int_{-1}^1 f(x, y, z) ~\text{d}x ~\text{d}y ~\text{d}z \approx \sum_{i=1}^N \sum_{j=1}^N \sum_{k=1}^N w_i\,w_j\,w_k \, f(x_i, y_i, z_i)
+```
+
+MeshIntegrals.jl uses [FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) to efficiently compute Gauss-Legendre quadrature nodes and weights at the time a `GaussLegendre` rule is constructed. Once constructed, these rules can be re-used to calculate multiple integrals to improve performance. Additionally, MeshIntegrals.jl uses an allocation-free summation routine that further improves performance by avoiding storing intermediate results.
+```julia
+rule = GaussLegendre(N)
+```
+
+By contrast to adaptive integration rules where compute times can sometimes vary significantly, this technique has the advantage that compute times are much more predictable. For example, calculating a one-dimensional integral whose integrand $f$ can be computed in 10 microseconds using a `GaussLegendre(100)` can be expected to take approximately 1 millisecond. However, this lacks many of the guard-rails present in adaptive routines: results are not automatically checked to ensure convergence, so care must be taken to ensure that an appropriate rule and order are chosen.
+
+## [H-Adaptive Cubature](@id hadaptivecubature)
+
+The `HAdaptiveCubature` type is used to specify an h-adaptive cubature rule, as implemented by [HCubature.jl](https://github.com/JuliaMath/HCubature.jl).
+
+```julia
+rule = HAdaptiveCubature()
+```
+
+All standard `HCubature.hcubature` keyword-argument options are supported. These can be specified when constructing the rule, where the `kwargs` in `HAdaptiveCubature(kwargs...)` is equivalent to `hcubature(f, a, b; kwargs...)`, e.g.:
+```julia
+rule = HAdaptiveCubature(order = 5, rtol = 1e-4)
+```