diff --git a/Project.toml b/Project.toml
index 7e9cfe3b..5a1d4711 100644
--- a/Project.toml
+++ b/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"
diff --git a/README.md b/README.md
index 68cc0ce0..1e51d48c 100644
--- a/README.md
+++ b/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
 
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 2ec3cf17..d830f578 100644
--- a/benchmark/benchmarks.jl
+++ b/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
 
diff --git a/docs/make.jl b/docs/make.jl
index 988736cb..fb2f4507 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -4,8 +4,7 @@ using MeshIntegrals
 # Dynamically set all files in subdirectories of the source directory to include all files in these subdirectories.
 # This way they don't need to be listed explicitly.
 SPECIALIZATIONS_FILES = joinpath.(Ref("specializations"),
-    readdir(joinpath(dirname(@__DIR__), "src",
-        "specializations")))
+    readdir(joinpath(dirname(@__DIR__), "src", "specializations")))
 
 makedocs(
     sitename = "MeshIntegrals.jl",
@@ -15,8 +14,9 @@ makedocs(
             "Support Matrix" => "supportmatrix.md",
             "Example Usage" => "usage.md"
         ],
-        "Derivations" => [
-            "Integrating a Triangle" => "triangle.md"
+        "Developer Notes" => [
+            "How it Works" => "how_it_works.md",
+            "Specializations" => "specializations.md"
         ],
         "Public API" => "api.md"
     ]
diff --git a/docs/src/how_it_works.md b/docs/src/how_it_works.md
new file mode 100644
index 00000000..3f230427
--- /dev/null
+++ b/docs/src/how_it_works.md
@@ -0,0 +1,76 @@
+# How it Works (By Example)
+
+## Example Problem
+
+Let $f$ be a function of position $\bar{r}$ in some space.
+```julia
+function f(r̄::Meshes.Point)
+    x, y, z = to(r̄)
+    ...
+end
+```
+
+Let the integration domain be the space (a ball) enclosed by a sphere centered on the origin with a radius of 5 meters.
+```julia
+center = Meshes.Point(0u"m", 0u"m", 0u"m")
+radius = 5.0u"m"
+ball = Meshes.Ball(center, radius)
+```
+
+This integral is often expressed abstractly as simply the following, where the triple integral signs and $\text{d}V$ indicate that the integration domain is some three-dimensional volume.
+```math
+\iiint f(\bar{r}) ~ \text{d}V
+```
+
+Integrals like this are often solved manually by selecting an appropriate coordinate system and limits that neatly represent the integration domain, e.g.
+```math
+\int_0^{\pi} \int_0^{2\pi} \int_0^{5} f(\bar{r}) ~ \text{d}\rho~\text{d}\theta~\text{d}\phi
+```
+
+This works great for simple geometries, but requires integration code that is geometry-specific. This package leverages parametric functions defined in Meshes.jl and differential forms to define integral methods that are general solutions for all geometries.
+
+## [Parametric Functions](@id how-parametric)
+
+Every supported `Meshes.Geometry` type is defined as having a parametric function that maps from a local parametric coordinate system to every point on the geometry. Curve-like geometries will have a single parametric dimension, surfaces will have two dimensions, and volumes will have three dimensions; this can be checked for a particular geometry via `Meshes.paramdim(geometry)`.
+
+For consistency across geometry types, with [some notable exceptions](@ref specializations), these parametric functions are defined to take coordinates inside a normalized range $[0,1]$. In the example case of `ball`, Meshes.jl defines a parametric function mapped in normalized spherical coordinates $(t_\rho, ~t_\theta, ~t_\phi)$. We find, then:
+```julia
+Meshes.paramdim(ball) == 3    # a volume
+
+ball(tρ, tθ, tφ)    # for args in range [0, 1], maps to a corresponding Meshes.Point
+
+ball(0, tθ, tφ) == center
+```
+
+In effect, we can now use the geometry itself as a function that maps from three normalized ($0 \le t \le 1$) arguments to every point on the geometry. For the sake of generalization, let this parametric function be called $g$.
+```math
+\text{g}: (t_1,~t_2,~t_3) ~\mapsto~ \text{Point}\big[ x, ~y, ~z \big]  
+```
+
+## 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(r̄) ~ \text{d}V = \iiint f(\bar{r}) ~ \bar{\text{d}t_1} \wedge \bar{\text{d}t_2} \wedge \bar{\text{d}t_3}
+```
+
+This resultant differential (volume) element is formed at each point in the integration domain by taking [the Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of the parametric function.
+```math
+\mathbf{J}_f = \begin{bmatrix} \bar{\text{d}t_1} & \bar{\text{d}t_2} & \bar{\text{d}t_3} \end{bmatrix}
+```
+where
+```math
+\bar{\text{d}t_n} = \frac{\partial}{\partial t_n} ~ \text{g}(t_1,~t_2,~t_3)
+```
+
+Each of these partial derivatives is a vector representing the direction that changing each parametric function argument will move the resultant point. The differential element ($E$) size is then calculated using geometric algebra as the magnitude of the exterior product ($\wedge$) of these three vectors.
+```math
+E(t_1,~t_2,~t_3) = \left\| \bar{\text{d}t_1} \wedge \bar{\text{d}t_2} \wedge \bar{\text{d}t_3} \right\|
+```
+
+Finally, we use the parametric function itself, $g$, as a map to all points $\bar{r}$ in the integration domain. Since `Meshes.Geometry` parametric functions all operate on normalized domains, we can now solve any volume integral as simply
+```math
+\int_0^1 \int_0^1 \int_0^1 f\Big(\text{g}\big(t_1,~t_2,~t_3\big)\Big) ~ E(t_1,~t_2,~t_3) ~ \text{d}t_1 ~ \text{d}t_2 ~ \text{d}t_3
+```
+
+This form of integral can be trivially generalized to support $n$-dimensional geometries in a form that enables the use of a wide range of numerical integration libraries.
diff --git a/docs/src/index.md b/docs/src/index.md
index 52edc831..f947eb64 100644
--- a/docs/src/index.md
+++ b/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)
diff --git a/docs/src/triangle.md b/docs/src/specializations.md
similarity index 61%
rename from docs/src/triangle.md
rename to docs/src/specializations.md
index a0e53d0c..dc873601 100644
--- a/docs/src/triangle.md
+++ b/docs/src/specializations.md
@@ -1,4 +1,18 @@
-# Integrating a Triangle
+# [Specializations](@id specializations)
+
+There are several notable exceptions to how Meshes.jl defines [parametric functions](@ref how-parametric).
+- `Meshes.ConeSurface` is essentially a composite type and has a parametric function that only maps the conical portion of the geometry, so the `Meshes.Disk` base element has to be integrated separately.
+- `Meshes.CylinderSurface` is essentially a composite type and has a parametric function that only maps the cylindrical portion of the geometry, so the `Meshes.Disk` element has to be integrated separately.
+- `Meshes.FrustumSurface` is essentially a composite type and has a parametric function that only maps the cylindrical portion of the geometry, so the top and bottom `Meshes.Disk` elements have to be integrated separately.
+- `Meshes.Line` represents a line of infinite length that passes through two points, and it has a parametric function that is valid on the domain $(-\infty, \infty)$.
+- `Meshes.Plane` represents a plane of infinite extent, and it has a parametric function that is valid on the domain $(-\infty, \infty)^2$.
+- `Meshes.Ray` represents a line that begins at a point and extends in a particular direction with infinite length, and it has a parametric function that is valid on the domain $[0, \infty)$.
+- `Meshes.Ring` is a composite type that lacks a parametric function, but can be decomposed into `Meshes.Segment`s and then integrated by adding together the individual integrals.
+- `Meshes.Rope` is a composite type that lacks a parametric function, but can be decomposed into `Meshes.Segment`s and then integrated by adding together the individual integrals.
+- `Meshes.Triangle` has a parametric function that takes coordinates on a 2D barycentric coordinate system. So, for `(::Meshes.Triangle)(t1, t2)`, the coordinates must obey: $t_1, t_2 \in [0,1]$ where $t_1 + t_2 \le 1$.
+- `Meshes.Tetrahedron` has a parametric function that takes coordinates on a 3D barycentric coordinate system. So, for `(::Meshes.Tetrahedron)(t1, t2)`, the coordinates must obey: $t_1, t_2, t_3 \in [0,1]$ where $t_1 + t_2 + t_3 \le 1$.
+
+## Triangle
 
 For a specified `Meshes.Triangle` surface with area $A$, let $u$ and $v$ be Barycentric coordinates that span the surface.
 ```math
diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl
index 19228769..5fbac40f 100644
--- a/src/MeshIntegrals.jl
+++ b/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, jacobian
+
+include("utils.jl")
 
 include("integration_rules.jl")
 export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature
diff --git a/src/differentiation.jl b/src/differentiation.jl
index acb580f1..78bd7952 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -1,42 +1,97 @@
 ################################################################################
-#                     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
+
+# Default constructors
+FiniteDifference{T}() where {T <: AbstractFloat} = FiniteDifference{T}(T(1e-6))
+FiniteDifference() = FiniteDifference{Float64}()
+
+"""
+    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
+
+# Future Support:
+#   struct AutoEnzyme <: DifferentiationMethod end
+#   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::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}
+) where {G <: Geometry, T <: AbstractFloat}
+    return jacobian(geometry, ts, _default_method(G))
+end
+
 function jacobian(
         geometry::Geometry,
-        ts::V;
-        ε = 1e-6
-) where {V <: Union{AbstractVector, Tuple}}
+        ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
+        diff_method::FiniteDifference
+) where {T <: AbstractFloat}
     Dim = Meshes.paramdim(geometry)
     if Dim != length(ts)
         throw(ArgumentError("ts must have same number of dimensions as geometry."))
     end
 
-    T = eltype(ts)
-    ε = T(ε)
-
     # Get the partial derivative along the n'th axis via finite difference
     #   approximation, where ts is the current parametric position
-    function ∂ₙr(ts, n)
+    function ∂ₙr(ts, n, ε)
         # 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
-        if ts[n] < T(0.01)
+        if ts[n] < 0.01
             # Right
             return (geometry(right...) - geometry(ts...)) / ε
-        elseif T(0.99) < ts[n]
+        elseif 0.99 < ts[n]
             # Left
             return (geometry(ts...) - geometry(left...)) / ε
         else
@@ -45,37 +100,7 @@ function jacobian(
         end
     end
 
-    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,)
+    return ntuple(n -> ∂ₙr(ts, n, T(diff_method.ε)), Dim)
 end
 
 ################################################################################
@@ -83,17 +108,25 @@ end
 ################################################################################
 
 """
-    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::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
+        diff_method::DifferentiationMethod = _default_method(G)
+) where {G <: Geometry, T <: AbstractFloat}
     # 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
diff --git a/src/integral.jl b/src/integral.jl
index d4925db7..1a1c7878 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -3,7 +3,7 @@
 ################################################################################
 
 """
-    integral(f, geometry[, rule]; FP=Float64)
+    integral(f, geometry[, rule]; diff_method=_default_method(geometry), FP=Float64)
 
 Numerically integrate a given function `f(::Point)` over the domain defined by
 a `geometry` using a particular numerical integration `rule` with floating point
@@ -12,12 +12,13 @@ precision of type `FP`.
 # Arguments
 - `f`: an integrand function with a method `f(::Meshes.Point)`
 - `geometry`: some `Meshes.Geometry` that defines the integration domain
-- `rule`: optionally, the `IntegrationRule` used for integration (by default `GaussKronrod()` in 1D and `HAdaptiveCubature()` else)
-- `FP`: optionally, the floating point precision desired (`Float64` by default)
+- `rule`: optionally, the `IntegrationRule` used for integration (by default
+`GaussKronrod()` in 1D and `HAdaptiveCubature()` else)
 
-Note that reducing `FP` below `Float64` will incur some loss of precision. By
-contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision
-(at the expense of longer runtimes).
+# Keyword Arguments
+- `diff_method::DifferentiationMethod = _default_method(geometry)`: the method to
+use for calculating Jacobians that are used to calculate differential elements
+- `FP = Float64`: the floating point precision desired.
 """
 function integral end
 
@@ -59,8 +60,9 @@ function _integral(
         f,
         geometry,
         rule::GaussLegendre;
-        FP::Type{T} = Float64
-) where {T <: AbstractFloat}
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_method(geometry)
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
     # Get Gauss-Legendre nodes and weights for a region [-1,1]^N
@@ -73,7 +75,7 @@ function _integral(
 
     function integrand((weights, nodes))
         ts = t.(nodes)
-        prod(weights) * f(geometry(ts...)) * differential(geometry, ts)
+        prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method)
     end
 
     return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes))
@@ -84,11 +86,12 @@ function _integral(
         f,
         geometry,
         rule::HAdaptiveCubature;
-        FP::Type{T} = Float64
-) where {T <: AbstractFloat}
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_method(geometry)
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
-    integrand(ts) = f(geometry(ts...)) * differential(geometry, ts)
+    integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, diff_method)
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
@@ -111,9 +114,10 @@ function _integral_gk_1d(
         f,
         geometry,
         rule::GaussKronrod;
-        FP::Type{T} = Float64
-) where {T <: AbstractFloat}
-    integrand(t) = f(geometry(t)) * differential(geometry, (t,))
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_method(geometry)
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
+    integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
@@ -121,9 +125,10 @@ function _integral_gk_2d(
         f,
         geometry2d,
         rule::GaussKronrod;
-        FP::Type{T} = Float64
-) where {T <: AbstractFloat}
-    integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v))
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_method(geometry2d)
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
+    integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method)
     ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1]
     return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
 end
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index c40b3366..fc5905b3 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -9,31 +9,43 @@
 #   keyword argument and pass through the specified algorithm choice.
 ################################################################################
 
+################################################################################
+#                              integral
+################################################################################
 """
     integral(f, curve::BezierCurve, rule = GaussKronrod();
-             FP=Float64, alg=Meshes.Horner())
-
-Like [`integral`](@ref) but integrates along the domain defined a `curve`. By
-default this uses Horner's method to improve performance when parameterizing
-the `curve` at the expense of a small loss of precision. Additional accuracy
-can be obtained by specifying the use of DeCasteljau's algorithm instead with
-`alg=Meshes.DeCasteljau()` but can come at a steep cost in memory allocations,
-especially for curves with a large number of control points.
+             diff_method=Analytical(), FP=Float64, alg=Meshes.Horner())
+
+Like [`integral`](@ref) but integrates along the domain defined by `curve`.
+
+# Arguments
+- `f`: an integrand function with a method `f(::Meshes.Point)`
+- `curve`: a `Meshes.BezierCurve` that defines the integration domain
+- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
+
+# Keyword Arguments
+- `diff_method::DifferentiationMethod = Analytical()`: the method to use for
+calculating Jacobians that are used to calculate differential elements
+- `FP = Float64`: the floating point precision desired
+- `alg = Meshes.Horner()`:  the method to use for parameterizing `curve`. Alternatively,
+`alg=Meshes.DeCasteljau()` can be specified for increased accuracy, but comes with a
+steep performance cost, especially for curves with a large number of control points.
 """
 function integral(
         f::F,
         curve::Meshes.BezierCurve,
         rule::GaussLegendre;
+        diff_method::DM = _default_method(curve),
         FP::Type{T} = Float64,
         alg::Meshes.BezierEvalMethod = Meshes.Horner()
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
     # Change of variables: x [-1,1] ↦ t [0,1]
     t(x) = FP(1 // 2) * x + FP(1 // 2)
     point(x) = curve(t(x), alg)
-    integrand(x) = f(point(x)) * differential(curve, (t(x),))
+    integrand(x) = f(point(x)) * differential(curve, (t(x),), diff_method)
 
     # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length]
     return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
@@ -43,11 +55,12 @@ function integral(
         f::F,
         curve::Meshes.BezierCurve,
         rule::GaussKronrod;
+        diff_method::DM = _default_method(curve),
         FP::Type{T} = Float64,
         alg::Meshes.BezierEvalMethod = Meshes.Horner()
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
     point(t) = curve(t, alg)
-    integrand(t) = f(point(t)) * differential(curve, (t,))
+    integrand(t) = f(point(t)) * differential(curve, (t,), diff_method)
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
@@ -55,11 +68,12 @@ function integral(
         f::F,
         curve::Meshes.BezierCurve,
         rule::HAdaptiveCubature;
+        diff_method::DM = _default_method(curve),
         FP::Type{T} = Float64,
         alg::Meshes.BezierEvalMethod = Meshes.Horner()
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
     point(t) = curve(t, alg)
-    integrand(ts) = f(point(only(ts))) * differential(curve, ts)
+    integrand(ts) = f(point(only(ts))) * differential(curve, ts, diff_method)
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
@@ -73,3 +87,40 @@ function integral(
     # Reapply units
     return value .* integrandunits
 end
+
+################################################################################
+#                               jacobian
+################################################################################
+
+function jacobian(
+        bz::Meshes.BezierCurve,
+        ts::V,
+        diff_method::Analytical
+) 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
+
+_has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true
diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index fbb398c9..4d4c4c89 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -11,8 +11,11 @@ function integral(
         f::F,
         line::Meshes.Line,
         rule::GaussLegendre;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Line, diff_method)
+
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
@@ -24,8 +27,8 @@ function integral(
     t′(x) = (1 + x^2) / (1 - x^2)^2
 
     # Integrate f along the Line
-    domainunits = _units(line(0))
-    integrand(x) = f(line(t(x))) * t′(x) * domainunits
+    differential(line, x) = t′(x) * _units(line(0))
+    integrand(x) = f(line(t(x))) * differential(line, x)
     return sum(w .* integrand(x) for (w, x) in zip(ws, xs))
 end
 
@@ -33,8 +36,11 @@ function integral(
         f::F,
         line::Meshes.Line,
         rule::GaussKronrod;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Line, diff_method)
+
     # Normalize the Line s.t. line(t) is distance t from origin point
     line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
 
@@ -48,8 +54,11 @@ function integral(
         f::F,
         line::Meshes.Line,
         rule::HAdaptiveCubature;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Line, diff_method)
+
     # Normalize the Line s.t. line(t) is distance t from origin point
     line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
 
@@ -58,8 +67,8 @@ function integral(
     t′(x) = (1 + x^2) / (1 - x^2)^2
 
     # Integrate f along the Line
-    domainunits = _units(line(0))
-    integrand(x::AbstractVector) = f(line(t(x[1]))) * t′(x[1]) * domainunits
+    differential(line, x) = t′(x) * _units(line(0))
+    integrand(x::AbstractVector) = f(line(t(x[1]))) * differential(line, x[1])
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
@@ -68,8 +77,14 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, FP[-1], FP[1]; rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, (-one(FP),), (one(FP),); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
 end
+
+################################################################################
+#                               jacobian
+################################################################################
+
+_has_analytical(::Type{T}) where {T <: Meshes.Line} = true
diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index 70243e8a..1b34405e 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -12,8 +12,11 @@ function integral(
         f::F,
         plane::Meshes.Plane,
         rule::GaussLegendre;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Plane, diff_method)
+
     # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]²
     xs, ws = _gausslegendre(FP, rule.n)
     wws = Iterators.product(ws, ws)
@@ -40,8 +43,11 @@ function integral(
         f::F,
         plane::Meshes.Plane,
         rule::GaussKronrod;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Plane, diff_method)
+
     # Normalize the Plane's orthogonal vectors
     uu = Meshes.unormalize(plane.u)
     uv = Meshes.unormalize(plane.v)
@@ -58,8 +64,11 @@ function integral(
         f::F,
         plane::Meshes.Plane,
         rule::HAdaptiveCubature;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Plane, diff_method)
+
     # Normalize the Plane's orthogonal vectors
     uu = Meshes.unormalize(plane.u)
     uv = Meshes.unormalize(plane.v)
@@ -87,3 +96,9 @@ function integral(
     # Reapply units
     return value .* integrandunits
 end
+
+################################################################################
+#                               jacobian
+################################################################################
+
+_has_analytical(::Type{T}) where {T <: Meshes.Plane} = true
diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index fce24530..9008b4a8 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -12,8 +12,11 @@ function integral(
         f::F,
         ray::Meshes.Ray,
         rule::GaussLegendre;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Ray, diff_method)
+
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
@@ -38,8 +41,11 @@ function integral(
         f::F,
         ray::Meshes.Ray,
         rule::GaussKronrod;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Ray, diff_method)
+
     # Normalize the Ray s.t. ray(t) is distance t from origin point
     ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v))
 
@@ -53,8 +59,11 @@ function integral(
         f::F,
         ray::Meshes.Ray,
         rule::HAdaptiveCubature;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Ray, diff_method)
+
     # Normalize the Ray s.t. ray(t) is distance t from origin point
     ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v))
 
@@ -78,3 +87,9 @@ function integral(
     # Reapply units
     return value .* integrandunits
 end
+
+################################################################################
+#                               jacobian
+################################################################################
+
+_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true
diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl
index 51cadf56..b2490ee2 100644
--- a/src/specializations/Ring.jl
+++ b/src/specializations/Ring.jl
@@ -8,6 +8,24 @@
 #   constituent Segment's, integrated separately, and then summed.
 ################################################################################
 
+"""
+    integral(f, ring::Ring, rule = GaussKronrod();
+             diff_method=FiniteDifference(), FP=Float64)
+
+Like [`integral`](@ref) but integrates along the domain defined by `ring`. The
+specified integration `rule` is applied independently to each segment formed by
+consecutive points in the Ring.
+
+# Arguments
+- `f`: an integrand function with a method `f(::Meshes.Point)`
+- `ring`: a `Ring` that defines the integration domain
+- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
+
+# Keyword Arguments
+- `diff_method::DifferentiationMethod = FiniteDifference()`: the method to use for
+calculating Jacobians that are used to calculate differential elements
+- `FP = Float64`: the floating point precision desired
+"""
 function integral(
         f::F,
         ring::Meshes.Ring,
diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index c68e82cc..7fcefa7f 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -8,6 +8,24 @@
 #   integrated separately, and then summed.
 ################################################################################
 
+"""
+    integral(f, rope::Rope, rule = GaussKronrod();
+             diff_method=FiniteDifference(), FP=Float64)
+
+Like [`integral`](@ref) but integrates along the domain defined by `rope`. The
+specified integration `rule` is applied independently to each segment formed by
+consecutive points in the Rope.
+
+# Arguments
+- `f`: an integrand function with a method `f(::Meshes.Point)`
+- `rope`: a `Rope` that defines the integration domain
+- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
+
+# Keyword Arguments
+- `diff_method::DifferentiationMethod = FiniteDifference()`: the method to use for
+calculating Jacobians that are used to calculate differential elements
+- `FP = Float64`: the floating point precision desired
+"""
 function integral(
         f::F,
         rope::Meshes.Rope,
diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index e4aac959..9fb2af9a 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -13,8 +13,9 @@ function integral(
         f::F,
         tetrahedron::Meshes.Tetrahedron,
         rule::GaussLegendre;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
     _error_unsupported_combination("Tetrahedron", "GaussLegendre")
 end
 
@@ -22,8 +23,11 @@ function integral(
         f::F,
         tetrahedron::Meshes.Tetrahedron,
         rule::GaussKronrod;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Tetrahedron, diff_method)
+
     nil = zero(FP)
     ∫uvw(u, v, w) = f(tetrahedron(u, v, w))
     ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), nil, FP(1 - v - w); rule.kwargs...)[1]
@@ -38,7 +42,14 @@ function integral(
         f::F,
         tetrahedron::Meshes.Tetrahedron,
         rule::HAdaptiveCubature;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
     _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature")
 end
+
+################################################################################
+#                               jacobian
+################################################################################
+
+_has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index d0b418bb..c63e715b 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -9,20 +9,15 @@
 #   documentation.
 ################################################################################
 
-"""
-    integral(f, triangle::Triangle, ::GaussLegendre; FP=Float64)
-
-Like [`integral`](@ref) but integrates over the surface of a `triangle`
-by transforming the triangle into a polar-barycentric coordinate system and
-using a Gauss-Legendre quadrature rule along each barycentric dimension of the
-triangle.
-"""
 function integral(
         f::F,
-        triangle::Meshes.Ngon{3},
+        triangle::Meshes.Triangle,
         rule::GaussLegendre;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Triangle, diff_method)
+
     # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2
     xs, ws = _gausslegendre(FP, rule.n)
     wws = Iterators.product(ws, ws)
@@ -53,18 +48,15 @@ function integral(
     return FP(π / 4) * Meshes.area(triangle) .* sum(integrand, zip(wws, xxs))
 end
 
-"""
-    integral(f, triangle::Triangle, ::GaussKronrod; FP=Float64)
-
-Like [`integral`](@ref) but integrates over the surface of a `triangle` using nested
-Gauss-Kronrod quadrature rules along each barycentric dimension of the triangle.
-"""
 function integral(
         f::F,
-        triangle::Meshes.Ngon{3},
+        triangle::Meshes.Triangle,
         rule::GaussKronrod;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Triangle, diff_method)
+
     # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0)
     #   i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du
     ∫u(u) = QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1]
@@ -74,19 +66,15 @@ function integral(
     return 2 * Meshes.area(triangle) .* ∫
 end
 
-"""
-    integral(f, triangle::Triangle, ::GaussKronrod; FP=Float64)
-
-Like [`integral`](@ref) but integrates over the surface of a `triangle` by
-transforming the triangle into a polar-barycentric coordinate system and using
-an h-adaptive cubature rule.
-"""
 function integral(
         f::F,
-        triangle::Meshes.Ngon{3},
+        triangle::Meshes.Triangle,
         rule::HAdaptiveCubature;
+        diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, T <: AbstractFloat}
+) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+    _guarantee_analytical(Meshes.Triangle, diff_method)
+
     # Integrate the Barycentric triangle by transforming it into polar coordinates
     #   with a modified radius
     #     R = r ( sinφ + cosφ )
@@ -104,3 +92,9 @@ function integral(
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
     return 2 * Meshes.area(triangle) .* ∫
 end
+
+################################################################################
+#                               jacobian
+################################################################################
+
+_has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true
diff --git a/src/utils.jl b/src/utils.jl
index 14de6e6b..cafbfb67 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -14,6 +14,34 @@ function _error_unsupported_combination(geometry, rule)
     throw(ArgumentError(msg))
 end
 
+################################################################################
+#                           DifferentiationMethod
+################################################################################
+
+# Throw an ArgumentError if Analytical() jacobian not defined for this type
+function _guarantee_analytical(
+        ::Type{G},
+        ::DifferentiationMethod
+) where {G <: Geometry}
+    throw(ArgumentError("$G geometries require kwarg diff_method = Analytical()"))
+end
+
+_guarantee_analytical(::Type{G}, ::Analytical) where {G <: Geometry} = nothing
+
+# Return whether a geometry type has jacobian methods defined
+_has_analytical(::Type{G}) where {G <: Geometry} = false
+_has_analytical(g::G) where {G <: Geometry} = _has_analytical(G)
+
+# Return the default DifferentiationMethod instance for a particular geometry type
+function _default_method(
+        g::Type{G}
+) where {G <: Geometry}
+    return _has_analytical(G) ? Analytical() : FiniteDifference()
+end
+
+# Return the default DifferentiationMethod instance for a particular geometry instance
+_default_method(g::G) where {G <: Geometry} = _default_method(G)
+
 ################################################################################
 #                        CliffordNumbers and Units
 ################################################################################
diff --git a/test/utils.jl b/test/utils.jl
index 45453a50..2a2d8a68 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -9,3 +9,35 @@
     p = Point(1.0u"cm", 2.0u"mm", 3.0u"m")
     @test MeshIntegrals._units(p) == u"m"
 end
+
+@testitem "DifferentiationMethod" setup=[Setup] begin
+    using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical
+
+    # _has_analytical of instances
+    bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)])
+    @test _has_analytical(bezier) == true
+    sphere = Sphere(Point(0, 0, 0), 1.0)
+    @test _has_analytical(sphere) == false
+
+    # _has_analytical of types
+    @test _has_analytical(Meshes.BezierCurve) == true
+    @test _has_analytical(Meshes.Line) == true
+    @test _has_analytical(Meshes.Plane) == true
+    @test _has_analytical(Meshes.Ray) == true
+    @test _has_analytical(Meshes.Sphere) == false
+    @test _has_analytical(Meshes.Tetrahedron) == true
+    @test _has_analytical(Meshes.Triangle) == true
+
+    # _guarantee_analytical
+    @test _guarantee_analytical(Meshes.Line, Analytical()) === nothing
+    @test_throws "Analytical" _guarantee_analytical(Meshes.Line, FiniteDifference())
+
+    # _default_method
+    @test _default_method(Meshes.BezierCurve) isa Analytical
+    @test _default_method(bezier) isa Analytical
+    @test _default_method(Meshes.Sphere) isa FiniteDifference
+    @test _default_method(sphere) isa FiniteDifference
+
+    # FiniteDifference
+    @test FiniteDifference().ε ≈ 1e-6
+end