add AD (Enzyme) support via EnzymeExt

.gitignore

@@ -22,3 +22,7 @@ docs/site/
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+
+# Development related
+.vscode
+dev/

Project.toml

@@ -1,6 +1,6 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
-version = "0.15.2"
+version = "0.16.0-DEV"
 
 [deps]
 CliffordNumbers = "3998ac73-6bd4-4031-8035-f167dd3ed523"
@@ -12,9 +12,16 @@ Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
+[weakdeps]
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+
+[extensions]
+MeshIntegralsEnzymeExt = "Enzyme"
+
 [compat]
 CliffordNumbers = "0.1.4"
 CoordRefSystems = "0.12, 0.13, 0.14, 0.15"
+Enzyme = "0.13.12"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 LinearAlgebra = "1"

README.md

@@ -20,27 +20,31 @@ These solvers have support for integrand functions that produce scalars, vectors
 
 ## Usage
 
+### Basic
+
 ```julia
 integral(f, geometry)
 ```
 Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry`. A default integration method will be automatically selected according to the geometry: `GaussKronrod()` for 1D, and `HAdaptiveCubature()` for all others.
 
 ```julia
-integral(f, geometry, rule, FP=Float64)
+integral(f, geometry, rule)
 ```
 Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`.
 
-Optionally, a fourth argument can be provided to specify the floating point precision level desired. This setting can be manipulated if your integrand function produces outputs with alternate floating point precision (e.g. `Float16`, `BigFloat`, etc) AND you'd prefer to avoid implicit type promotions.
+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.
+
+### Aliases
 
 ```julia
 lineintegral(f, geometry)
 surfaceintegral(f, geometry)
 volumeintegral(f, geometry)
 ```
-Alias functions are provided for convenience. These are simply wrappers for `integral` that first validate that the provided `geometry` has the expected number of parametric/manifold dimensions. Like with `integral` in the examples above, the `rule` can also be specified as a third-argument.
-- `lineintegral` for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc)
-- `surfaceintegral` for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
-- `volumeintegral` for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)
+Alias functions are provided for convenience. These are simply wrappers for `integral` that also validate that the provided `geometry` has the expected number of parametric dimensions. Like with `integral`, a `rule` can also optionally be specified as a third argument.
+- `lineintegral` is used for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc)
+- `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
+- `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)
 
 # Demo
 

benchmark/benchmarks.jl

@@ -28,7 +28,7 @@ geometries = (
 SUITE["Integrals"] = let s = BenchmarkGroup()
     for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
         n1, n2, N = geometry.name, "$(int.name) $(rule.name)", rule.N
-        s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule) evals=N
+        s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
     end
     s
 end
@@ -41,8 +41,8 @@ sphere = Sphere(Point(0, 0, 0), 1.0)
 differential = MeshIntegrals.differential
 
 SUITE["Differentials"] = let s = BenchmarkGroup()
-    s["Jacobian"] = @benchmarkable jacobian($sphere, $(0.5, 0.5)) evals=1000
-    s["Differential"] = @benchmarkable differential($sphere, $(0.5, 0.5)) evals=1000
+    s["Jacobian"] = @benchmarkable jacobian($sphere, $(0.5, 0.5)) evals=10
+    s["Differential"] = @benchmarkable differential($sphere, $(0.5, 0.5)) evals=10
     s
 end
 

docs/make.jl

@@ -12,6 +12,7 @@ makedocs(
     pages = [
         "Home" => [
             "About" => "index.md",
+            "How it Works" => "how_it_works.md",
             "Support Matrix" => "supportmatrix.md",
             "Example Usage" => "usage.md"
         ],

docs/src/how_it_works.md

@@ -0,0 +1,43 @@
+# How it Works (Work in Progress)
+
+This page will explain how this package works by example...
+
+Let $f$ be a function to be integrated throughout the volume bounded by a unit sphere. This integral is often expressed as simply
+```math
+\iiint f(x, y, z) ~ \text{d}V
+```
+
+## Parametric Functions
+
+**!TODO! update segment to Ball**
+
+Every supported `Meshes.Geometry` type is defined as having a parametric function that maps from a local coordinate system to every point on the geometry. For example,
+```julia
+a = Meshes.Point(0, 0)
+b = Meshes.Point(2, 4)
+segment = Meshes.Segment(a, b)
+```
+defines a line segment beginning at point `a` and ending at point `b`. As a geometry with one parametric dimension (i.e. `paramdim(segment) == 1`), it can be treated as a function with one argument that returns a corresponding point on the segment.
+```julia
+segment(0) == a
+segment(0.5) == Point(1, 2)
+segment(1) == b
+```
+
+... where $t$ is a parametric coordinate used to generate points in the domain.
+
+
+
+## Differential Forms
+
+Using differential forms, the general solution for integrating a geometry with three parametric dimensions ($t_1$, $t_2$, and $t_3$) is
+```math
+\iiint f(x, y, z) ~ \left\| \text{d}t_1 \wedge \text{d}t_2 \wedge \text{d}t_3 \right\|
+```
+
+Since `Meshes.Geometry`s parametric functions have arguments that are defined on the domain $[0,1]$, this is equivalent to
+```math
+\int_0^1 \int_0^1 \int_0^1 f(x, y, z) ~ \left\| \text{d}t_1 \wedge \text{d}t_2 \wedge \text{d}t_3 \right\|
+```
+
+For every point in the integration domain where the integrand function is evaluated, the differential element $\text{d}t_1 \wedge \text{d}t_2$ is calculated using [the Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of the parametric function. For a two-dimensional surface this Jacobian consists of two vectors, each pointing in the direction that the parametric function's output point will move by changing each input argument. The differential element, then, is simply the magnitude of the exterior product of these vectors.

docs/src/index.md

@@ -10,5 +10,38 @@
 [![Coveralls](https://coveralls.io/repos/github/JuliaGeometry/MeshIntegrals.jl/badge.svg?branch=main)](https://coveralls.io/github/JuliaGeometry/MeshIntegrals.jl?branch=main)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
-**MeshIntegrals.jl** is a Julia library that leverages differential forms to implement fast
-and easy numerical integration of field equations over geometric domains.
+
+**MeshIntegrals.jl** uses differential forms to enable fast and easy numerical integration of arbitrary integrand functions over domains defined via [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) geometries. This is achieved using:
+- Gauss-Legendre quadrature rules from [**FastGaussQuadrature.jl**](https://github.com/JuliaApproximation/FastGaussQuadrature.jl): `GaussLegendre(n)`
+- H-adaptive Gauss-Kronrod quadrature rules from [**QuadGK.jl**](https://github.com/JuliaMath/QuadGK.jl): `GaussKronrod(kwargs...)`
+- H-adaptive cubature rules from [**HCubature.jl**](https://github.com/JuliaMath/HCubature.jl): `HAdaptiveCubature(kwargs...)`
+
+These solvers have support for integrand functions that produce scalars, vectors, and [**Unitful.jl**](https://github.com/PainterQubits/Unitful.jl) `Quantity` types. While HCubature.jl does not natively support `Quantity` type integrands, this package provides a compatibility layer to enable this feature.
+
+## Usage
+
+### Basic
+
+```julia
+integral(f, geometry)
+```
+Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry`. A default integration method will be automatically selected according to the geometry: `GaussKronrod()` for 1D, and `HAdaptiveCubature()` for all others.
+
+```julia
+integral(f, geometry, rule)
+```
+Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`.
+
+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.
+
+### Aliases
+
+```julia
+lineintegral(f, geometry)
+surfaceintegral(f, geometry)
+volumeintegral(f, geometry)
+```
+Alias functions are provided for convenience. These are simply wrappers for `integral` that also validate that the provided `geometry` has the expected number of parametric dimensions. Like with `integral`, a `rule` can also optionally be specified as a third argument.
+- `lineintegral` is used for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc)
+- `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
+- `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)

ext/MeshIntegralsEnzymeExt.jl

@@ -0,0 +1,15 @@
+module MeshIntegralsEnzymeExt
+
+using MeshIntegrals: MeshIntegrals, AutoEnzyme
+using Meshes: Meshes
+using Enzyme: Enzyme
+
+function MeshIntegrals.jacobian(
+        geometry::Meshes.Geometry,
+        ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
+        ::AutoEnzyme
+) where {T <: AbstractFloat}
+    return Meshes.to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...))
+end
+
+end

src/MeshIntegrals.jl

@@ -9,10 +9,10 @@ import LinearAlgebra
 import QuadGK
 import Unitful
 
-include("utils.jl")
-
 include("differentiation.jl")
-export jacobian
+export DifferentiationMethod, Analytical, FiniteDifference, AutoEnzyme, jacobian
+
+include("utils.jl")
 
 include("integration_rules.jl")
 export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature

src/differentiation.jl

@@ -1,35 +1,98 @@
 ################################################################################
-#                     Jacobian and Differential Elements
+#                          DifferentiationMethods
 ################################################################################
 
 """
-    jacobian(geometry, ts; ε=1e-6)
+    DifferentiationMethod
 
-Calculate the Jacobian of a geometry at some parametric point `ts` using a simple
-finite-difference approximation with step size `ε`.
+A category of types used to specify the desired method for calculating derivatives.
+Derivatives are used to form Jacobian matrices when calculating the differential
+element size throughout the integration region.
+
+See also [`FiniteDifference`](@ref), [`Analytical`](@ref).
+"""
+abstract type DifferentiationMethod end
+
+"""
+    FiniteDifference(ε=1e-6)
+
+Use to specify use of a finite-difference approximation method with a step size
+of `ε` for calculating derivatives.
+"""
+struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod
+    ε::T
+end
+
+# If ε not specified, default to 1e-6
+FiniteDifference() = FiniteDifference(1e-6)
+
+"""
+    Analytical()
+
+Use to specify use of analytically-derived solutions for calculating derivatives.
+These solutions are currently defined only for a subset of geometry types.
+
+# Supported Geometries:
+- `BezierCurve`
+- `Line`
+- `Plane`
+- `Ray`
+- `Tetrahedron`
+- `Triangle`
+"""
+struct Analytical <: DifferentiationMethod end
+
+"""
+    AutoEnzyme()
+
+Use to specify use of the Enzyme.jl for calculating derivatives.
+"""
+struct AutoEnzyme <: DifferentiationMethod end
+
+# Future Support:
+#   struct AutoZygote <: DifferentiationMethod end
+
+################################################################################
+#                                  Jacobian
+################################################################################
+
+"""
+    jacobian(geometry, ts[, diff_method])
+
+Calculate the Jacobian of a `geometry`'s parametric function at some point `ts`.
+Optionally, direct the use of a particular differentiation method `diff_method`; by default
+use analytic solutions where possible and finite difference approximations
+otherwise.
 
 # Arguments
 - `geometry`: some `Meshes.Geometry` of N parametric dimensions
 - `ts`: a parametric point specified as a vector or tuple of length N
-- `ε`: step size to use for the finite-difference approximation
+- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use
 """
+function jacobian(
+        geometry::G,
+        ts::V
+) where {G <: Geometry, V <: Union{AbstractVector, Tuple}}
+    return jacobian(geometry, ts, _default_method(G))
+end
+
 function jacobian(
         geometry::Geometry,
-        ts::V;
-        ε = 1e-6
+        ts::V,
+        diff_method::FiniteDifference
 ) where {V <: Union{AbstractVector, Tuple}}
     Dim = Meshes.paramdim(geometry)
     if Dim != length(ts)
         throw(ArgumentError("ts must have same number of dimensions as geometry."))
     end
 
     T = eltype(ts)
-    ε = T(ε)
+    ε = T(diff_method.ε)
 
     # Get the partial derivative along the n'th axis via finite difference
     #   approximation, where ts is the current parametric position
     function ∂ₙr(ts, n)
-        # Build left/right parametric coordinates with non-allocating iterators 
+        # Build left/right parametric coordinates with non-allocating iterators
         left = Iterators.map(((i, t),) -> i == n ? t - ε : t, enumerate(ts))
         right = Iterators.map(((i, t),) -> i == n ? t + ε : t, enumerate(ts))
         # Select orientation of finite-diff
@@ -48,52 +111,30 @@ function jacobian(
     return ntuple(n -> ∂ₙr(ts, n), Dim)
 end
 
-function jacobian(
-        bz::Meshes.BezierCurve,
-        ts::V
-) where {V <: Union{AbstractVector, Tuple}}
-    t = only(ts)
-    # Parameter t restricted to domain [0,1] by definition
-    if t < 0 || t > 1
-        throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
-    end
-
-    # Aliases
-    P = bz.controls
-    N = Meshes.degree(bz)
-
-    # Ensure that this implementation is tractible: limited by ability to calculate
-    #   binomial(N, N/2) without overflow. It's possible to extend this range by
-    #   converting N to a BigInt, but this results in always returning BigFloat's.
-    N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.")
-
-    # Generator for Bernstein polynomial functions
-    B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i)
-
-    # Derivative = N Σ_{i=0}^{N-1} sigma(i)
-    #   P indices adjusted for Julia 1-based array indexing
-    sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1])
-    derivative = N .* sum(sigma, 0:(N - 1))
-
-    return (derivative,)
-end
-
 ################################################################################
 #                          Differential Elements
 ################################################################################
 
 """
-    differential(geometry, ts)
+    differential(geometry, ts[, diff_method])
 
 Calculate the differential element (length, area, volume, etc) of the parametric
-function for `geometry` at arguments `ts`.
+function for `geometry` at arguments `ts`. Optionally, direct the use of a
+particular differentiation method `diff_method`; by default use analytic solutions where
+possible and finite difference approximations otherwise.
+
+# Arguments
+- `geometry`: some `Meshes.Geometry` of N parametric dimensions
+- `ts`: a parametric point specified as a vector or tuple of length N
+- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use
 """
 function differential(
-        geometry::Geometry,
-        ts::V
-) where {V <: Union{AbstractVector, Tuple}}
+        geometry::G,
+        ts::V,
+        diff_method::DifferentiationMethod = _default_method(G)
+) where {G <: Geometry, V <: Union{AbstractVector, Tuple}}
     # Calculate the Jacobian, convert Vec -> KVector
-    J = jacobian(geometry, ts)
+    J = jacobian(geometry, ts, diff_method)
     J_kvecs = Iterators.map(_kvector, J)
 
     # Extract units from Geometry type