From f1107cf8dd446eb19fa964bde25c31b8d5a0e495 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Thu, 12 Dec 2024 22:17:11 -0500
Subject: [PATCH 01/24] add Enzyme as a potential differentiation method for
 the jacobian

---
 .gitignore                             |  4 +++
 Project.toml                           |  7 +++++
 ext/MeshIntegralsEnzymeExt.jl          | 15 ++++++++++
 src/MeshIntegrals.jl                   |  2 +-
 src/differentiation.jl                 | 33 +++++++++++++++++----
 src/integral.jl                        | 10 +++++--
 src/specializations/BezierCurve.jl     |  7 +++--
 src/specializations/CylinderSurface.jl | 11 ++++---
 src/utils.jl                           | 11 ++++++-
 test/Project.toml                      |  1 +
 test/combinations.jl                   | 41 ++++++++++++++++----------
 11 files changed, 111 insertions(+), 31 deletions(-)
 create mode 100644 ext/MeshIntegralsEnzymeExt.jl

diff --git a/.gitignore b/.gitignore
index 29126e47..d4dad176 100644
--- a/.gitignore
+++ b/.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
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index dc58d615..cfc86403 100644
--- a/Project.toml
+++ b/Project.toml
@@ -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.9"
 CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
+Enzyme = "0.13.22"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 LinearAlgebra = "1"
diff --git a/ext/MeshIntegralsEnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl
new file mode 100644
index 00000000..ccc6f36d
--- /dev/null
+++ b/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
diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl
index cbc7260f..82da7fc4 100644
--- a/src/MeshIntegrals.jl
+++ b/src/MeshIntegrals.jl
@@ -10,7 +10,7 @@ import QuadGK
 import Unitful
 
 include("differentiation.jl")
-export DifferentiationMethod, FiniteDifference, jacobian
+export DifferentiationMethod, FiniteDifference, AutoEnzyme, jacobian
 
 include("utils.jl")
 
diff --git a/src/differentiation.jl b/src/differentiation.jl
index 4433fd49..a05b4062 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -9,7 +9,7 @@ A category of types used to specify the desired method for calculating derivativ
 Derivatives are used to form Jacobian matrices when calculating the differential
 element size throughout the integration region.
 
-See also [`FiniteDifference`](@ref).
+See also [`FiniteDifference`](@ref), [`Analytical`](@ref).
 """
 abstract type DifferentiationMethod end
 
@@ -23,12 +23,33 @@ struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod
     ε::T
 end
 
-# Default constructors
-FiniteDifference{T}() where {T <: AbstractFloat} = FiniteDifference{T}(T(1e-6))
-FiniteDifference() = FiniteDifference{Float64}()
+# 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 AutoEnzyme <: DifferentiationMethod end
 #   struct AutoZygote <: DifferentiationMethod end
 
 ################################################################################
@@ -68,7 +89,7 @@ function jacobian(
     # 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
diff --git a/src/integral.jl b/src/integral.jl
index 34d9695e..41c27c0f 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -3,7 +3,7 @@
 ################################################################################
 
 """
-    integral(f, geometry[, rule]; diff_method=_default_method(geometry), FP=Float64)
+    integral(f, geometry[, rule]; diff_method=_default_diff_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
@@ -16,7 +16,7 @@ precision of type `FP`.
 `GaussKronrod()` in 1D and `HAdaptiveCubature()` else)
 
 # Keyword Arguments
-- `diff_method::DifferentiationMethod = _default_method(geometry)`: the method to
+- `diff_method::DifferentiationMethod = _default_diff_method(geometry)`: the method to
 use for calculating Jacobians that are used to calculate differential elements
 - `FP = Float64`: the floating point precision desired.
 """
@@ -44,6 +44,8 @@ function _integral(
         FP::Type{T} = Float64,
         diff_method::DM = _default_diff_method(geometry)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
+    _check_diff_method_support(geometry, diff_method)
+
     # Implementation depends on number of parametric dimensions over which to integrate
     N = Meshes.paramdim(geometry)
     if N == 1
@@ -72,6 +74,8 @@ function _integral(
         FP::Type{T} = Float64,
         diff_method::DM = _default_diff_method(geometry)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
+    _check_diff_method_support(geometry, diff_method)
+
     N = Meshes.paramdim(geometry)
 
     # Get Gauss-Legendre nodes and weights of type FP for a region [-1,1]ᴺ
@@ -101,6 +105,8 @@ function _integral(
         FP::Type{T} = Float64,
         diff_method::DM = _default_diff_method(geometry)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
+    _check_diff_method_support(geometry, diff_method)
+
     N = Meshes.paramdim(geometry)
 
     integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, diff_method)
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index cdae2ec6..6007cb67 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -36,13 +36,16 @@ function integral(
         curve::Meshes.BezierCurve,
         rule::IntegrationRule;
         alg::Meshes.BezierEvalMethod = Meshes.Horner(),
+        diff_method::DM = _default_diff_method(curve),
         kwargs...
-)
+) where {DM <: DifferentiationMethod}
+    _check_diff_method_support(curve, diff_method)
+
     # Generate a _ParametricGeometry whose parametric function auto-applies the alg kwarg
     param_curve = _ParametricGeometry(_parametric(curve, alg), Meshes.paramdim(curve))
 
     # Integrate the _ParametricGeometry using the standard methods
-    return _integral(f, param_curve, rule; kwargs...)
+    return _integral(f, param_curve, rule; diff_method = diff_method, kwargs...)
 end
 
 ################################################################################
diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl
index 3e79d0be..627b4018 100644
--- a/src/specializations/CylinderSurface.jl
+++ b/src/specializations/CylinderSurface.jl
@@ -12,18 +12,21 @@ function integral(
         f,
         cyl::Meshes.CylinderSurface,
         rule::I;
+        diff_method::DM = _default_diff_method(cyl),
         kwargs...
-) where {I <: IntegrationRule}
+) where {I <: IntegrationRule, DM <: DifferentiationMethod}
+    _check_diff_method_support(cyl, diff_method)
+
     # The generic method only parametrizes the sides
-    sides = _integral(f, cyl, rule; kwargs...)
+    sides = _integral(f, cyl, rule; diff_method = diff_method, kwargs...)
 
     # Integrate the Disk at the top
     disk_top = Meshes.Disk(cyl.top, cyl.radius)
-    top = _integral(f, disk_top, rule; kwargs...)
+    top = _integral(f, disk_top, rule; diff_method = diff_method, kwargs...)
 
     # Integrate the Disk at the bottom
     disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
-    bottom = _integral(f, disk_bottom, rule; kwargs...)
+    bottom = _integral(f, disk_bottom, rule; diff_method = diff_method, kwargs...)
 
     return sides + top + bottom
 end
diff --git a/src/utils.jl b/src/utils.jl
index 6cf5ee14..1f6b4864 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -12,6 +12,8 @@ end
 #                           DifferentiationMethod
 ################################################################################
 
+_check_diff_method_support(::Geometry, ::DifferentiationMethod) = nothing
+
 # Return the default DifferentiationMethod instance for a particular geometry type
 function _default_diff_method(
         g::Type{G}
@@ -20,7 +22,14 @@ function _default_diff_method(
 end
 
 # Return the default DifferentiationMethod instance for a particular geometry instance
-_default_diff_method(g::G) where {G <: Geometry} = _default_diff_method(G)
+_default_diff_method(::G) where {G <: Geometry} = _default_diff_method(G)
+
+non_enzyme_types = (:BezierCurve, :CylinderSurface, :Cylinder, :ParametrizedCurve)
+for geometry_type in non_enzyme_types
+    @eval function _check_diff_method_support(::Meshes.$geometry_type, ::AutoEnzyme)
+        throw(ArgumentError("Differentiation method AutoEnzyme not supported for $(string(Meshes.$geometry_type))."))
+    end
+end
 
 ################################################################################
 #                           Numerical Tools
diff --git a/test/Project.toml b/test/Project.toml
index e38225db..10f2a671 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
diff --git a/test/combinations.jl b/test/combinations.jl
index b53420bd..cb2b2593 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -18,6 +18,7 @@ This file includes tests for:
     using Meshes
     using MeshIntegrals
     using Unitful
+    using Enzyme
 
     # Used for testing callable objects as integrand functions
     struct Callable{F <: Function}
@@ -43,27 +44,27 @@ This file includes tests for:
         gausslegendre::Bool
         hadaptivecubature::Bool
         # DifferentiationMethods
-        # autoenzyme::Bool
+        autoenzyme::Bool
     end
 
     # Shortcut constructor for geometries with typical support structure
-    function SupportStatus(geometry::Geometry)
+    function SupportStatus(geometry::Geometry, autoenzyme = true)
         if paramdim(geometry) == 1
             aliases = Bool.((1, 0, 0))
             rules = Bool.((1, 1, 1))
-            return SupportStatus(aliases..., rules...)
+            return SupportStatus(aliases..., rules..., autoenzyme)
         elseif paramdim(geometry) == 2
             aliases = Bool.((0, 1, 0))
             rules = Bool.((1, 1, 1))
-            return SupportStatus(aliases..., rules...)
+            return SupportStatus(aliases..., rules..., autoenzyme)
         elseif paramdim(geometry) == 3
             aliases = Bool.((0, 0, 1))
             rules = Bool.((0, 1, 1))
-            return SupportStatus(aliases..., rules...)
+            return SupportStatus(aliases..., rules..., autoenzyme)
         else
             aliases = Bool.((0, 0, 0))
             rules = Bool.((0, 1, 1))
-            return SupportStatus(aliases..., rules...)
+            return SupportStatus(aliases..., rules..., autoenzyme)
         end
     end
 
@@ -110,15 +111,19 @@ This file includes tests for:
             end
         end # for
 
-        #=
         iter_diff_methods = (
             (supports.autoenzyme, AutoEnzyme()),
         )
 
         for (supported, diff_method) in iter_diff_methods
-            @test integral(testable.integrand, testable.geometry; diff_method=diff_method)≈sol rtol=rtol
-        end
-        =#
+            if supported
+                @test integral(
+                    testable.integrand, testable.geometry; diff_method = diff_method)≈testable.solution rtol=rtol
+            else
+                @test_throws "not supported" integral(
+                    testable.integrand, testable.geometry; diff_method = diff_method)
+            end
+        end # for
     end # function
 end #testsnippet
 
@@ -180,7 +185,8 @@ end
 
     # Package and run tests
     testable = TestableGeometry(integrand, curve, solution)
-    runtests(testable; rtol = 0.5e-2)
+    supports = SupportStatus(curve, false)
+    runtests(testable, supports; rtol = 0.5e-2)
 end
 
 @testitem "Meshes.Box 1D" setup=[Combinations] begin
@@ -321,7 +327,8 @@ end
 
     # Package and run tests
     testable = TestableGeometry(integrand, cyl, solution)
-    runtests(testable)
+    supports = SupportStatus(cyl, false)
+    runtests(testable, supports)
 end
 
 @testitem "Meshes.CylinderSurface" setup=[Combinations] begin
@@ -336,7 +343,8 @@ end
 
     # Package and run tests
     testable = TestableGeometry(integrand, cyl, solution)
-    runtests(testable)
+    supports = SupportStatus(cyl, false)
+    runtests(testable, supports)
 end
 
 @testitem "Meshes.Disk" setup=[Combinations] begin
@@ -476,9 +484,12 @@ end
 
         # Package and run tests
         testable_cart = TestableGeometry(integrand, curve_cart, solution)
-        runtests(testable_cart)
+        supports_cart = SupportStatus(curve_cart, false)
+        runtests(testable_cart, supports_cart)
+
         testable_polar = TestableGeometry(integrand, curve_polar, solution)
-        runtests(testable_polar)
+        supports_polar = SupportStatus(curve_polar, false)
+        runtests(testable_polar, supports_polar)
     end
 end
 

From c776dc94a752254fb951a241c48754cbefbbd8fc Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 09:10:37 -0500
Subject: [PATCH 02/24] refactor check for enzyme support

---
 Project.toml           |  4 ++--
 src/differentiation.jl | 18 +-----------------
 src/utils.jl           | 21 +++++++++++++--------
 test/Project.toml      |  3 ++-
 test/combinations.jl   | 30 +++++++++++++++++++-----------
 test/utils.jl          |  4 ++--
 6 files changed, 39 insertions(+), 41 deletions(-)

diff --git a/Project.toml b/Project.toml
index cfc86403..6884be18 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,11 +21,11 @@ MeshIntegralsEnzymeExt = "Enzyme"
 [compat]
 CliffordNumbers = "0.1.9"
 CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
-Enzyme = "0.13.22"
+Enzyme = "0.13"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 LinearAlgebra = "1"
-Meshes = "0.50, 0.51, 0.52"
+Meshes = "0.51.20, 0.52"
 QuadGK = "2.1.1"
 Unitful = "1.19"
 julia = "1.9"
diff --git a/src/differentiation.jl b/src/differentiation.jl
index a05b4062..154124b2 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -9,7 +9,7 @@ A category of types used to specify the desired method for calculating derivativ
 Derivatives are used to form Jacobian matrices when calculating the differential
 element size throughout the integration region.
 
-See also [`FiniteDifference`](@ref), [`Analytical`](@ref).
+See also [`FiniteDifference`](@ref), [`AutoEnzyme`](@ref).
 """
 abstract type DifferentiationMethod end
 
@@ -26,22 +26,6 @@ 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()
 
diff --git a/src/utils.jl b/src/utils.jl
index 1f6b4864..f3d2f3e4 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -12,25 +12,30 @@ end
 #                           DifferentiationMethod
 ################################################################################
 
+supports_autoenzyme(::Type{<:Meshes.Geometry}) = true
+supports_autoenzyme(::Type{<:Meshes.BezierCurve}) = false
+supports_autoenzyme(::Type{<:Meshes.CylinderSurface}) = false
+supports_autoenzyme(::Type{<:Meshes.Cylinder}) = false
+supports_autoenzyme(::Type{<:Meshes.ParametrizedCurve}) = false
+supports_autoenzyme(::G) where {G <: Geometry} = supports_autoenzyme(G)
+
 _check_diff_method_support(::Geometry, ::DifferentiationMethod) = nothing
+function _check_diff_method_support(geometry::Geometry, ::AutoEnzyme)
+    if !supports_autoenzyme(geometry)
+        throw(ArgumentError("Differentiation method AutoEnzyme not supported for this geometry."))
+    end
+end
 
 # Return the default DifferentiationMethod instance for a particular geometry type
 function _default_diff_method(
         g::Type{G}
 ) where {G <: Geometry}
-    return FiniteDifference()
+    supports_autoenzyme(g) ? AutoEnzyme() : FiniteDifference()
 end
 
 # Return the default DifferentiationMethod instance for a particular geometry instance
 _default_diff_method(::G) where {G <: Geometry} = _default_diff_method(G)
 
-non_enzyme_types = (:BezierCurve, :CylinderSurface, :Cylinder, :ParametrizedCurve)
-for geometry_type in non_enzyme_types
-    @eval function _check_diff_method_support(::Meshes.$geometry_type, ::AutoEnzyme)
-        throw(ArgumentError("Differentiation method AutoEnzyme not supported for $(string(Meshes.$geometry_type))."))
-    end
-end
-
 ################################################################################
 #                           Numerical Tools
 ################################################################################
diff --git a/test/Project.toml b/test/Project.toml
index 10f2a671..6aeb90c4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -14,9 +14,10 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 Aqua = "0.7, 0.8"
 CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
+Enzyme = "0.13"
 ExplicitImports = "1.6.0"
 LinearAlgebra = "1"
-Meshes = "0.50, 0.51, 0.52"
+Meshes = "0.51.20, 0.52"
 SpecialFunctions = "2"
 TestItemRunner = "1"
 TestItems = "1"
diff --git a/test/combinations.jl b/test/combinations.jl
index cb2b2593..4a268d02 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -48,7 +48,8 @@ This file includes tests for:
     end
 
     # Shortcut constructor for geometries with typical support structure
-    function SupportStatus(geometry::Geometry, autoenzyme = true)
+    function SupportStatus(
+            geometry::Geometry, autoenzyme = MeshIntegrals.supports_autoenzyme(geometry))
         if paramdim(geometry) == 1
             aliases = Bool.((1, 0, 0))
             rules = Bool.((1, 1, 1))
@@ -183,10 +184,12 @@ end
     end
     solution = 2π * u"Ω*m"
 
+    # Enzyme support
+    @test MeshIntegrals.supports_autoenzyme(curve) == false
+
     # Package and run tests
     testable = TestableGeometry(integrand, curve, solution)
-    supports = SupportStatus(curve, false)
-    runtests(testable, supports; rtol = 0.5e-2)
+    runtests(testable; rtol = 0.5e-2)
 end
 
 @testitem "Meshes.Box 1D" setup=[Combinations] begin
@@ -325,10 +328,12 @@ end
     integrand(p) = 1.0u"A"
     solution = Meshes.measure(cyl) * u"A"
 
+    # Enzyme support
+    @test MeshIntegrals.supports_autoenzyme(cyl) == false
+
     # Package and run tests
     testable = TestableGeometry(integrand, cyl, solution)
-    supports = SupportStatus(cyl, false)
-    runtests(testable, supports)
+    runtests(testable)
 end
 
 @testitem "Meshes.CylinderSurface" setup=[Combinations] begin
@@ -341,10 +346,12 @@ end
     integrand(p) = 1.0u"A"
     solution = Meshes.measure(cyl) * u"A"
 
+    # Enzyme support
+    @test MeshIntegrals.supports_autoenzyme(cyl) == false
+
     # Package and run tests
     testable = TestableGeometry(integrand, cyl, solution)
-    supports = SupportStatus(cyl, false)
-    runtests(testable, supports)
+    runtests(testable)
 end
 
 @testitem "Meshes.Disk" setup=[Combinations] begin
@@ -482,14 +489,15 @@ end
         end
         solution = 2π * radius * exp(-radius^2) * u"A*m"
 
+        # Enzyme support
+        @test MeshIntegrals.supports_autoenzyme(curve_cart) == false
+
         # Package and run tests
         testable_cart = TestableGeometry(integrand, curve_cart, solution)
-        supports_cart = SupportStatus(curve_cart, false)
-        runtests(testable_cart, supports_cart)
+        runtests(testable_cart)
 
         testable_polar = TestableGeometry(integrand, curve_polar, solution)
-        supports_polar = SupportStatus(curve_polar, false)
-        runtests(testable_polar, supports_polar)
+        runtests(testable_polar)
     end
 end
 
diff --git a/test/utils.jl b/test/utils.jl
index 947688a7..9f3e254b 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -26,8 +26,8 @@ end
 @testitem "Differentiation" setup=[Utils] begin
     # _default_diff_method
     sphere = Sphere(Point(0, 0, 0), 1.0)
-    @test _default_diff_method(Meshes.Sphere) isa FiniteDifference
-    @test _default_diff_method(sphere) isa FiniteDifference
+    @test _default_diff_method(Meshes.Sphere) isa AutoEnzyme
+    @test _default_diff_method(sphere) isa AutoEnzyme
 
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6

From 9f3aadc244fa9dc131712e0695459415d7f9057a Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 09:39:03 -0500
Subject: [PATCH 03/24] add FP to _default_diff_method

---
 ext/MeshIntegralsEnzymeExt.jl          |  4 ++++
 src/differentiation.jl                 |  4 ++--
 src/integral.jl                        | 10 +++++-----
 src/specializations/BezierCurve.jl     |  7 ++++---
 src/specializations/CylinderSurface.jl | 11 ++++++-----
 src/utils.jl                           | 14 ++++++++++----
 test/floating_point_types.jl           |  1 +
 test/utils.jl                          |  5 +++--
 8 files changed, 35 insertions(+), 21 deletions(-)

diff --git a/ext/MeshIntegralsEnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl
index ccc6f36d..339b8b64 100644
--- a/ext/MeshIntegralsEnzymeExt.jl
+++ b/ext/MeshIntegralsEnzymeExt.jl
@@ -9,6 +9,10 @@ function MeshIntegrals.jacobian(
         ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
         ::AutoEnzyme
 ) where {T <: AbstractFloat}
+    Dim = Meshes.paramdim(geometry)
+    if Dim != length(ts)
+        throw(ArgumentError("ts must have same number of dimensions as geometry."))
+    end
     return Meshes.to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...))
 end
 
diff --git a/src/differentiation.jl b/src/differentiation.jl
index 154124b2..5c63b8e9 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -57,7 +57,7 @@ function jacobian(
         geometry::G,
         ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}
 ) where {G <: Geometry, T <: AbstractFloat}
-    return jacobian(geometry, ts, _default_diff_method(G))
+    return jacobian(geometry, ts, _default_diff_method(G, T))
 end
 
 function jacobian(
@@ -112,7 +112,7 @@ possible and finite difference approximations otherwise.
 function differential(
         geometry::G,
         ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
-        diff_method::DifferentiationMethod = _default_diff_method(G)
+        diff_method::DifferentiationMethod = _default_diff_method(G, T)
 ) where {G <: Geometry, T <: AbstractFloat}
     J = Iterators.map(_KVector, jacobian(geometry, ts, diff_method))
     return LinearAlgebra.norm(foldl(∧, J))
diff --git a/src/integral.jl b/src/integral.jl
index 41c27c0f..1365aa93 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -3,7 +3,7 @@
 ################################################################################
 
 """
-    integral(f, geometry[, rule]; diff_method=_default_diff_method(geometry), FP=Float64)
+    integral(f, geometry[, rule]; diff_method=_default_diff_method(geometry, FP), 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
@@ -16,7 +16,7 @@ precision of type `FP`.
 `GaussKronrod()` in 1D and `HAdaptiveCubature()` else)
 
 # Keyword Arguments
-- `diff_method::DifferentiationMethod = _default_diff_method(geometry)`: the method to
+- `diff_method::DifferentiationMethod = _default_diff_method(geometry, FP)`: the method to
 use for calculating Jacobians that are used to calculate differential elements
 - `FP = Float64`: the floating point precision desired.
 """
@@ -42,7 +42,7 @@ function _integral(
         geometry,
         rule::GaussKronrod;
         FP::Type{T} = Float64,
-        diff_method::DM = _default_diff_method(geometry)
+        diff_method::DM = _default_diff_method(geometry, FP)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _check_diff_method_support(geometry, diff_method)
 
@@ -72,7 +72,7 @@ function _integral(
         geometry,
         rule::GaussLegendre;
         FP::Type{T} = Float64,
-        diff_method::DM = _default_diff_method(geometry)
+        diff_method::DM = _default_diff_method(geometry, FP)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _check_diff_method_support(geometry, diff_method)
 
@@ -103,7 +103,7 @@ function _integral(
         geometry,
         rule::HAdaptiveCubature;
         FP::Type{T} = Float64,
-        diff_method::DM = _default_diff_method(geometry)
+        diff_method::DM = _default_diff_method(geometry, FP)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _check_diff_method_support(geometry, diff_method)
 
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index 6007cb67..5a8778d5 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -36,16 +36,17 @@ function integral(
         curve::Meshes.BezierCurve,
         rule::IntegrationRule;
         alg::Meshes.BezierEvalMethod = Meshes.Horner(),
-        diff_method::DM = _default_diff_method(curve),
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_diff_method(curve, FP),
         kwargs...
-) where {DM <: DifferentiationMethod}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _check_diff_method_support(curve, diff_method)
 
     # Generate a _ParametricGeometry whose parametric function auto-applies the alg kwarg
     param_curve = _ParametricGeometry(_parametric(curve, alg), Meshes.paramdim(curve))
 
     # Integrate the _ParametricGeometry using the standard methods
-    return _integral(f, param_curve, rule; diff_method = diff_method, kwargs...)
+    return _integral(f, param_curve, rule; diff_method = diff_method, FP = FP, kwargs...)
 end
 
 ################################################################################
diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl
index 627b4018..c359e3af 100644
--- a/src/specializations/CylinderSurface.jl
+++ b/src/specializations/CylinderSurface.jl
@@ -12,21 +12,22 @@ function integral(
         f,
         cyl::Meshes.CylinderSurface,
         rule::I;
-        diff_method::DM = _default_diff_method(cyl),
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_diff_method(cyl, FP),
         kwargs...
-) where {I <: IntegrationRule, DM <: DifferentiationMethod}
+) where {I <: IntegrationRule, DM <: DifferentiationMethod, T <: AbstractFloat}
     _check_diff_method_support(cyl, diff_method)
 
     # The generic method only parametrizes the sides
-    sides = _integral(f, cyl, rule; diff_method = diff_method, kwargs...)
+    sides = _integral(f, cyl, rule; diff_method = diff_method, FP = FP, kwargs...)
 
     # Integrate the Disk at the top
     disk_top = Meshes.Disk(cyl.top, cyl.radius)
-    top = _integral(f, disk_top, rule; diff_method = diff_method, kwargs...)
+    top = _integral(f, disk_top, rule; diff_method = diff_method, FP = FP, kwargs...)
 
     # Integrate the Disk at the bottom
     disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
-    bottom = _integral(f, disk_bottom, rule; diff_method = diff_method, kwargs...)
+    bottom = _integral(f, disk_bottom, rule; diff_method = diff_method, FP = FP, kwargs...)
 
     return sides + top + bottom
 end
diff --git a/src/utils.jl b/src/utils.jl
index f3d2f3e4..6556743e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -28,13 +28,19 @@ end
 
 # Return the default DifferentiationMethod instance for a particular geometry type
 function _default_diff_method(
-        g::Type{G}
-) where {G <: Geometry}
-    supports_autoenzyme(g) ? AutoEnzyme() : FiniteDifference()
+        g::Type{G}, FP::Type{T}
+) where {G <: Geometry, T <: AbstractFloat}
+    if supports_autoenzyme(g) && FP <: Union{Float32, Float64}
+        AutoEnzyme()
+    else
+        FiniteDifference()
+    end
 end
 
 # Return the default DifferentiationMethod instance for a particular geometry instance
-_default_diff_method(::G) where {G <: Geometry} = _default_diff_method(G)
+function _default_diff_method(::G, ::Type{T}) where {G <: Geometry, T <: AbstractFloat}
+    _default_diff_method(G, T)
+end
 
 ################################################################################
 #                           Numerical Tools
diff --git a/test/floating_point_types.jl b/test/floating_point_types.jl
index 7af87bd5..d1f6dbe5 100644
--- a/test/floating_point_types.jl
+++ b/test/floating_point_types.jl
@@ -7,6 +7,7 @@
     using LinearAlgebra: norm
     using Meshes
     using Unitful
+    using Enzyme
 
     baseatol = Dict(
         Float32 => 0.01f0,
diff --git a/test/utils.jl b/test/utils.jl
index 9f3e254b..038fcd81 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -26,8 +26,9 @@ end
 @testitem "Differentiation" setup=[Utils] begin
     # _default_diff_method
     sphere = Sphere(Point(0, 0, 0), 1.0)
-    @test _default_diff_method(Meshes.Sphere) isa AutoEnzyme
-    @test _default_diff_method(sphere) isa AutoEnzyme
+    @test _default_diff_method(Meshes.Sphere, Float64) isa AutoEnzyme
+    @test _default_diff_method(sphere, Float64) isa AutoEnzyme
+    @test _default_diff_method(sphere, BigFloat) isa FiniteDifference
 
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6

From dfe12a1c6692b8b5a0f6fe9c56e134603927ac76 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:04:04 -0500
Subject: [PATCH 04/24] add `using Enzyme` to benchmarks.jl

---
 benchmark/benchmarks.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 10ec4757..3535e840 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -3,6 +3,7 @@ using LinearAlgebra
 using Meshes
 using MeshIntegrals
 using Unitful
+using Enzyme
 
 const SUITE = BenchmarkGroup()
 

From 539500b20696e4d6e530f2d4dec31cfbb6933378 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:08:14 -0500
Subject: [PATCH 05/24] update CoordRefSystems.jl compat

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 6884be18..a4659137 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,7 +20,7 @@ MeshIntegralsEnzymeExt = "Enzyme"
 
 [compat]
 CliffordNumbers = "0.1.9"
-CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
+CoordRefSystems = "0.15, 0.16"
 Enzyme = "0.13"
 FastGaussQuadrature = "1"
 HCubature = "1.5"

From 875356fa74bbfa2e0b92975f5cf3a0a44a87940e Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:10:06 -0500
Subject: [PATCH 06/24] add Enzyme to Benchmark Project.toml

---
 benchmark/Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 96f687ac..3ac2c6cc 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -1,11 +1,13 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 BenchmarkTools = "1.5"
+Enzyme = "0.13"
 LinearAlgebra = "1"
 Meshes = "0.50, 0.51, 0.52"
 Unitful = "1.19"

From c22113797b3fcb0fd63cfc4cecce84748ee902d4 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:41:22 -0500
Subject: [PATCH 07/24] fix Meshes compat in Benchmark Project.toml

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 benchmark/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 3ac2c6cc..58968dc8 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -9,6 +9,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 BenchmarkTools = "1.5"
 Enzyme = "0.13"
 LinearAlgebra = "1"
-Meshes = "0.50, 0.51, 0.52"
+Meshes = "0.51.20, 0.52
 Unitful = "1.19"
 julia = "1.9"

From 8ecff85b55c2c5f53230b04789e76b7eec2f9bac Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:41:55 -0500
Subject: [PATCH 08/24] use import Enzyme, not using Enzyme

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 benchmark/benchmarks.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 3535e840..14896a03 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -3,7 +3,7 @@ using LinearAlgebra
 using Meshes
 using MeshIntegrals
 using Unitful
-using Enzyme
+import Enzyme
 
 const SUITE = BenchmarkGroup()
 

From 76999ac0a2b0c72a77b0f368262f87e96f1a4bdb Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:45:21 -0500
Subject: [PATCH 09/24] fix typo in Benchmarks Project.toml

---
 benchmark/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 58968dc8..709e5f12 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -9,6 +9,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 BenchmarkTools = "1.5"
 Enzyme = "0.13"
 LinearAlgebra = "1"
-Meshes = "0.51.20, 0.52
+Meshes = "0.51.20, 0.52"
 Unitful = "1.19"
 julia = "1.9"

From bd8174dde5cc2ef01da276d7cb639ef896f6d873 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:53:13 -0500
Subject: [PATCH 10/24] remove Meshes version check in combinations.jl

---
 test/combinations.jl | 50 ++++++++++++++++++++------------------------
 1 file changed, 23 insertions(+), 27 deletions(-)

diff --git a/test/combinations.jl b/test/combinations.jl
index 4a268d02..e6a8c295 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -469,36 +469,32 @@ end
 end
 
 @testitem "ParametrizedCurve" setup=[Combinations] begin
-    # ParametrizedCurve has been added in Meshes v0.51.20
-    # If the version is specified as minimal compat bound in the Project.toml, the downgrade test fails
-    if pkgversion(Meshes) >= v"0.51.20"
-        using CoordRefSystems: Polar
-
-        # Geometries
-        # Parametrize a circle centered on origin with specified radius
-        radius = 4.4
-        curve_cart = ParametrizedCurve(
-            t -> Point(radius * cos(t), radius * sin(t)), (0.0, 2π))
-        curve_polar = ParametrizedCurve(t -> Point(Polar(radius, t)), (0.0, 2π))
-
-        # Integrand & Solution
-        function integrand(p::Meshes.Point)
-            ur = norm(to(p))
-            r = ustrip(u"m", ur)
-            exp(-r^2) * u"A"
-        end
-        solution = 2π * radius * exp(-radius^2) * u"A*m"
-
-        # Enzyme support
-        @test MeshIntegrals.supports_autoenzyme(curve_cart) == false
+    using CoordRefSystems: Polar
 
-        # Package and run tests
-        testable_cart = TestableGeometry(integrand, curve_cart, solution)
-        runtests(testable_cart)
+    # Geometries
+    # Parametrize a circle centered on origin with specified radius
+    radius = 4.4
+    curve_cart = ParametrizedCurve(
+        t -> Point(radius * cos(t), radius * sin(t)), (0.0, 2π))
+    curve_polar = ParametrizedCurve(t -> Point(Polar(radius, t)), (0.0, 2π))
 
-        testable_polar = TestableGeometry(integrand, curve_polar, solution)
-        runtests(testable_polar)
+    # Integrand & Solution
+    function integrand(p::Meshes.Point)
+        ur = norm(to(p))
+        r = ustrip(u"m", ur)
+        exp(-r^2) * u"A"
     end
+    solution = 2π * radius * exp(-radius^2) * u"A*m"
+
+    # Enzyme support
+    @test MeshIntegrals.supports_autoenzyme(curve_cart) == false
+
+    # Package and run tests
+    testable_cart = TestableGeometry(integrand, curve_cart, solution)
+    runtests(testable_cart)
+    
+    testable_polar = TestableGeometry(integrand, curve_polar, solution)
+    runtests(testable_polar)
 end
 
 @testitem "Meshes.Plane" setup=[Combinations] begin

From 1786576c9d2dc34e86de236f41eeee0164fdc1a9 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 10:57:05 -0500
Subject: [PATCH 11/24] Apply format suggestion

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 test/combinations.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/combinations.jl b/test/combinations.jl
index e6a8c295..7faaca0d 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -492,7 +492,6 @@ end
     # Package and run tests
     testable_cart = TestableGeometry(integrand, curve_cart, solution)
     runtests(testable_cart)
-    
     testable_polar = TestableGeometry(integrand, curve_polar, solution)
     runtests(testable_polar)
 end

From 374715441706cf6cc05830eebc8648d8982580e1 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 11:42:24 -0500
Subject: [PATCH 12/24] Update test/Project.toml

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index 6aeb90c4..c5287142 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,7 +13,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Aqua = "0.7, 0.8"
-CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
+CoordRefSystems = "0.15, 0.16"
 Enzyme = "0.13"
 ExplicitImports = "1.6.0"
 LinearAlgebra = "1"

From 0d669a0a4983a3c94a6da9affb3f103eed5c5738 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Date: Fri, 13 Dec 2024 19:19:29 +0100
Subject: [PATCH 13/24] Bump compat of Enzyme to v0.13.19

---
 Project.toml           | 2 +-
 benchmark/Project.toml | 2 +-
 test/Project.toml      | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index a4659137..1e9cbb14 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,7 +21,7 @@ MeshIntegralsEnzymeExt = "Enzyme"
 [compat]
 CliffordNumbers = "0.1.9"
 CoordRefSystems = "0.15, 0.16"
-Enzyme = "0.13"
+Enzyme = "0.13.19"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 LinearAlgebra = "1"
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 709e5f12..6eec58b7 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -7,7 +7,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 BenchmarkTools = "1.5"
-Enzyme = "0.13"
+Enzyme = "0.13.19"
 LinearAlgebra = "1"
 Meshes = "0.51.20, 0.52"
 Unitful = "1.19"
diff --git a/test/Project.toml b/test/Project.toml
index c5287142..0861ae09 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -14,7 +14,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 Aqua = "0.7, 0.8"
 CoordRefSystems = "0.15, 0.16"
-Enzyme = "0.13"
+Enzyme = "0.13.19"
 ExplicitImports = "1.6.0"
 LinearAlgebra = "1"
 Meshes = "0.51.20, 0.52"

From f74aa98b4566a5222b6f68c47956347d87b65e92 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 13 Dec 2024 17:35:06 -0500
Subject: [PATCH 14/24] test supports_autoenzyme to combinations; test both
 backends for wrong dims

---
 test/combinations.jl | 14 ++------------
 test/utils.jl        |  4 +++-
 2 files changed, 5 insertions(+), 13 deletions(-)

diff --git a/test/combinations.jl b/test/combinations.jl
index 7faaca0d..4885c2d1 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -120,9 +120,11 @@ This file includes tests for:
             if supported
                 @test integral(
                     testable.integrand, testable.geometry; diff_method = diff_method)≈testable.solution rtol=rtol
+                @test MeshIntegrals.supports_autoenzyme(testable.geometry) == true
             else
                 @test_throws "not supported" integral(
                     testable.integrand, testable.geometry; diff_method = diff_method)
+                @test MeshIntegrals.supports_autoenzyme(testable.geometry) == false
             end
         end # for
     end # function
@@ -184,9 +186,6 @@ end
     end
     solution = 2π * u"Ω*m"
 
-    # Enzyme support
-    @test MeshIntegrals.supports_autoenzyme(curve) == false
-
     # Package and run tests
     testable = TestableGeometry(integrand, curve, solution)
     runtests(testable; rtol = 0.5e-2)
@@ -328,9 +327,6 @@ end
     integrand(p) = 1.0u"A"
     solution = Meshes.measure(cyl) * u"A"
 
-    # Enzyme support
-    @test MeshIntegrals.supports_autoenzyme(cyl) == false
-
     # Package and run tests
     testable = TestableGeometry(integrand, cyl, solution)
     runtests(testable)
@@ -346,9 +342,6 @@ end
     integrand(p) = 1.0u"A"
     solution = Meshes.measure(cyl) * u"A"
 
-    # Enzyme support
-    @test MeshIntegrals.supports_autoenzyme(cyl) == false
-
     # Package and run tests
     testable = TestableGeometry(integrand, cyl, solution)
     runtests(testable)
@@ -486,9 +479,6 @@ end
     end
     solution = 2π * radius * exp(-radius^2) * u"A*m"
 
-    # Enzyme support
-    @test MeshIntegrals.supports_autoenzyme(curve_cart) == false
-
     # Package and run tests
     testable_cart = TestableGeometry(integrand, curve_cart, solution)
     runtests(testable_cart)
diff --git a/test/utils.jl b/test/utils.jl
index 038fcd81..e2fd6be0 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -3,6 +3,7 @@
     using Meshes
     using MeshIntegrals: _default_diff_method, _parametric, _units, _zeros, _ones
     using Unitful
+    using Enzyme
 end
 
 @testitem "Utilities" setup=[Utils] begin
@@ -35,7 +36,8 @@ end
 
     # Test jacobian with wrong number of parametric coordinates
     box = Box(Point(0, 0), Point(1, 1))
-    @test_throws ArgumentError jacobian(box, zeros(3))
+    @test_throws ArgumentError MeshIntegrals.jacobian(box, zeros(3), FiniteDifference())
+    @test_throws ArgumentError MeshIntegrals.jacobian(box, zeros(3), AutoEnzyme())
 end
 
 @testitem "_ParametricGeometry" setup=[Utils] begin

From 08b9bf2aca339619457c2688791f81f0c764aa3d Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Fri, 13 Dec 2024 23:18:03 -0500
Subject: [PATCH 15/24] Restore recently-updated FiniteDifference constructors

---
 src/differentiation.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/differentiation.jl b/src/differentiation.jl
index 5c63b8e9..58425ad8 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -23,8 +23,9 @@ struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod
     ε::T
 end
 
-# If ε not specified, default to 1e-6
-FiniteDifference() = FiniteDifference(1e-6)
+# Default constructors
+FiniteDifference{T}() where {T <: AbstractFloat} = FiniteDifference{T}(T(1e-6))
+FiniteDifference() = FiniteDifference{Float64}()
 
 """
     AutoEnzyme()

From 45949a127afb5fc21204c5190097ac1c45ef292c Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Fri, 13 Dec 2024 23:33:45 -0500
Subject: [PATCH 16/24] Add docstrings, formatting

---
 src/utils.jl | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index 6556743e..84ebfdc7 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -12,6 +12,12 @@ end
 #                           DifferentiationMethod
 ################################################################################
 
+"""
+    supports_autoenzyme(geometry)
+
+Return whether a geometry (or geometry type) has a parametric function that can be
+differentiated with Enzyme. See GitHub Issue #154 for more information.
+"""
 supports_autoenzyme(::Type{<:Meshes.Geometry}) = true
 supports_autoenzyme(::Type{<:Meshes.BezierCurve}) = false
 supports_autoenzyme(::Type{<:Meshes.CylinderSurface}) = false
@@ -19,14 +25,24 @@ supports_autoenzyme(::Type{<:Meshes.Cylinder}) = false
 supports_autoenzyme(::Type{<:Meshes.ParametrizedCurve}) = false
 supports_autoenzyme(::G) where {G <: Geometry} = supports_autoenzyme(G)
 
+"""
+    _check_diff_method_support(::Geometry, ::DifferentiationMethod) -> nothing
+
+Throw an error if incompatible geometry-diff_method combination detected.
+"""
 _check_diff_method_support(::Geometry, ::DifferentiationMethod) = nothing
 function _check_diff_method_support(geometry::Geometry, ::AutoEnzyme)
     if !supports_autoenzyme(geometry)
-        throw(ArgumentError("Differentiation method AutoEnzyme not supported for this geometry."))
+        throw(ArgumentError("AutoEnzyme not supported for this geometry."))
     end
 end
 
-# Return the default DifferentiationMethod instance for a particular geometry type
+"""
+    _default_diff_method(geometry, FP)
+
+Return an instance of the default DifferentiationMethod for a particular geometry
+(or geometry type) and floating point type.
+"""
 function _default_diff_method(
         g::Type{G}, FP::Type{T}
 ) where {G <: Geometry, T <: AbstractFloat}
@@ -37,7 +53,6 @@ function _default_diff_method(
     end
 end
 
-# Return the default DifferentiationMethod instance for a particular geometry instance
 function _default_diff_method(::G, ::Type{T}) where {G <: Geometry, T <: AbstractFloat}
     _default_diff_method(G, T)
 end

From e05c8adf501f1f3cb97389bb3f0865a688e9a0eb Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Fri, 13 Dec 2024 23:37:40 -0500
Subject: [PATCH 17/24] Formatting

---
 test/combinations.jl | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/test/combinations.jl b/test/combinations.jl
index 4885c2d1..6becd2bd 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -48,21 +48,25 @@ This file includes tests for:
     end
 
     # Shortcut constructor for geometries with typical support structure
-    function SupportStatus(
-            geometry::Geometry, autoenzyme = MeshIntegrals.supports_autoenzyme(geometry))
-        if paramdim(geometry) == 1
+    function SupportStatus(g::Geometry, autoenzyme = MeshIntegrals.supports_autoenzyme(g))
+        N = Meshes.paramdim(g)
+        if N == 1
+            # line/curve
             aliases = Bool.((1, 0, 0))
             rules = Bool.((1, 1, 1))
             return SupportStatus(aliases..., rules..., autoenzyme)
-        elseif paramdim(geometry) == 2
+        elseif N == 2
+            # surface
             aliases = Bool.((0, 1, 0))
             rules = Bool.((1, 1, 1))
             return SupportStatus(aliases..., rules..., autoenzyme)
-        elseif paramdim(geometry) == 3
+        elseif N == 3
+            # volume
             aliases = Bool.((0, 0, 1))
             rules = Bool.((0, 1, 1))
             return SupportStatus(aliases..., rules..., autoenzyme)
         else
+            # ≥4D
             aliases = Bool.((0, 0, 0))
             rules = Bool.((0, 1, 1))
             return SupportStatus(aliases..., rules..., autoenzyme)

From 8d194138e7534764d3c222283923eb59af9b7727 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Sat, 14 Dec 2024 00:24:14 -0500
Subject: [PATCH 18/24] Add test for two-arg jacobian

---
 test/utils.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/utils.jl b/test/utils.jl
index e2fd6be0..08518a3c 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -34,6 +34,10 @@ end
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6
 
+    # Two-argument jacobian
+    segment = Segment(Point(0), Point(1))
+    @test jacobian(segment, (0.5,)) == Vec(1)
+
     # Test jacobian with wrong number of parametric coordinates
     box = Box(Point(0, 0), Point(1, 1))
     @test_throws ArgumentError MeshIntegrals.jacobian(box, zeros(3), FiniteDifference())

From 570a53e070500f9332c0ff67fc880d3648ee98ef Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Sat, 14 Dec 2024 00:39:18 -0500
Subject: [PATCH 19/24] Use rest of MeshIntegrals namespace

---
 test/utils.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/utils.jl b/test/utils.jl
index 08518a3c..2c0dae02 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -1,6 +1,7 @@
 @testsnippet Utils begin
     using LinearAlgebra: norm
     using Meshes
+    using MeshIntegrals
     using MeshIntegrals: _default_diff_method, _parametric, _units, _zeros, _ones
     using Unitful
     using Enzyme

From 71215f13fe17a5bf6e452f4bc7ba5c862b5d283c Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Sat, 14 Dec 2024 00:40:56 -0500
Subject: [PATCH 20/24] Disambiguate use of jacobian

---
 test/utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/utils.jl b/test/utils.jl
index 2c0dae02..3b9ca9c0 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -37,7 +37,7 @@ end
 
     # Two-argument jacobian
     segment = Segment(Point(0), Point(1))
-    @test jacobian(segment, (0.5,)) == Vec(1)
+    @test MeshIntegrals.jacobian(segment, (0.5,)) == Vec(1)
 
     # Test jacobian with wrong number of parametric coordinates
     box = Box(Point(0, 0), Point(1, 1))

From d2eff8319b5a0a7bd8e11e81683e1c9b59c59729 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Date: Sat, 14 Dec 2024 14:01:21 +0100
Subject: [PATCH 21/24] fix test

---
 test/utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/utils.jl b/test/utils.jl
index 3b9ca9c0..9f570a97 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -37,7 +37,7 @@ end
 
     # Two-argument jacobian
     segment = Segment(Point(0), Point(1))
-    @test MeshIntegrals.jacobian(segment, (0.5,)) == Vec(1)
+    @test MeshIntegrals.jacobian(segment, (0.5,)) == (Vec(1),)
 
     # Test jacobian with wrong number of parametric coordinates
     box = Box(Point(0, 0), Point(1, 1))

From c4f6b9280fbae1f6baf80e3fb8afd775b31a6dc2 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sat, 14 Dec 2024 10:07:09 -0500
Subject: [PATCH 22/24] use `import Enzyme`

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 test/combinations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/combinations.jl b/test/combinations.jl
index 6becd2bd..d7db64fb 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -18,7 +18,7 @@ This file includes tests for:
     using Meshes
     using MeshIntegrals
     using Unitful
-    using Enzyme
+    import Enzyme
 
     # Used for testing callable objects as integrand functions
     struct Callable{F <: Function}

From 7f25b8a9f96891f644de5f7618f0252c1fa31f0d Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sat, 14 Dec 2024 10:07:18 -0500
Subject: [PATCH 23/24] use `import Enzyme`

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 test/utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/utils.jl b/test/utils.jl
index 9f570a97..2ce8a6ab 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -4,7 +4,7 @@
     using MeshIntegrals
     using MeshIntegrals: _default_diff_method, _parametric, _units, _zeros, _ones
     using Unitful
-    using Enzyme
+    import Enzyme
 end
 
 @testitem "Utilities" setup=[Utils] begin

From 65dfe9c94c979daeb24033b8109ab5837c26eb23 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sat, 14 Dec 2024 10:07:42 -0500
Subject: [PATCH 24/24] remove unneeded MeshIntegrals.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 test/utils.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/utils.jl b/test/utils.jl
index 2ce8a6ab..5db3d779 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -41,8 +41,8 @@ end
 
     # Test jacobian with wrong number of parametric coordinates
     box = Box(Point(0, 0), Point(1, 1))
-    @test_throws ArgumentError MeshIntegrals.jacobian(box, zeros(3), FiniteDifference())
-    @test_throws ArgumentError MeshIntegrals.jacobian(box, zeros(3), AutoEnzyme())
+    @test_throws ArgumentError jacobian(box, zeros(3), FiniteDifference())
+    @test_throws ArgumentError jacobian(box, zeros(3), AutoEnzyme())
 end
 
 @testitem "_ParametricGeometry" setup=[Utils] begin