From 533dbfa3608e87fe0f4a8e9b4866bab190ee0914 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 11:38:29 -0500
Subject: [PATCH 01/22] Remove Function requirement

---
 src/integral.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/integral.jl b/src/integral.jl
index de69db5d..e750220f 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -24,7 +24,7 @@ function integral end
 
 # If only f and geometry are specified, select default rule
 function integral(
-        f::Function,
+        f,
         geometry::Geometry,
         rule::I = Meshes.paramdim(geometry) == 1 ? GaussKronrod() : HAdaptiveCubature();
         kwargs...

From 65b1547dee810583a6e3e99355c73b58679f6428 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 11:40:15 -0500
Subject: [PATCH 02/22] Remove Function requirement

---
 src/integral_aliases.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/integral_aliases.jl b/src/integral_aliases.jl
index 0981fec4..88e01bf3 100644
--- a/src/integral_aliases.jl
+++ b/src/integral_aliases.jl
@@ -15,7 +15,7 @@ Rule types available:
 - [`HAdaptiveCubature`](@ref)
 """
 function lineintegral(
-        f::Function,
+        f,
         geometry::Geometry,
         rule::IntegrationRule = GaussKronrod();
         kwargs...
@@ -47,7 +47,7 @@ Algorithm types available:
 - [`HAdaptiveCubature`](@ref) (default)
 """
 function surfaceintegral(
-        f::Function,
+        f,
         geometry::Geometry,
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
@@ -79,7 +79,7 @@ Algorithm types available:
 - [`HAdaptiveCubature`](@ref) (default)
 """
 function volumeintegral(
-        f::Function,
+        f,
         geometry::Geometry,
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...

From 8534d4d47db82210cf2203d4446a0b9a3ed44304 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 11:52:55 -0500
Subject: [PATCH 03/22] Remove Function requirement

---
 src/specializations/BezierCurve.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index c3540cbc..bd465978 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -32,13 +32,13 @@ calculating Jacobians that are used to calculate differential elements
 steep performance cost, especially for curves with a large number of control points.
 """
 function integral(
-        f::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, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
@@ -52,26 +52,26 @@ function integral(
 end
 
 function integral(
-        f::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, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     point(t) = curve(t, alg)
     integrand(t) = f(point(t)) * differential(curve, (t,), diff_method)
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
 function integral(
-        f::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, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     point(t) = curve(t, alg)
     integrand(ts) = f(point(only(ts))) * differential(curve, ts, diff_method)
 

From df7bbf484adb1df0507009bd1fac4f5ae991531c Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 11:53:30 -0500
Subject: [PATCH 04/22] Remove Function requirement

---
 src/specializations/ConeSurface.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/specializations/ConeSurface.jl b/src/specializations/ConeSurface.jl
index d22aaa77..7d0071b1 100644
--- a/src/specializations/ConeSurface.jl
+++ b/src/specializations/ConeSurface.jl
@@ -9,11 +9,11 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         cone::Meshes.ConeSurface,
         rule::I;
         kwargs...
-) where {F <: Function, I <: IntegrationRule}
+) where {I <: IntegrationRule}
     # The generic method only parameterizes the sides
     sides = _integral(f, cone, rule; kwargs...)
 

From 93c5398a403787be2e7cd735f5b4e713804a95a9 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 11:54:03 -0500
Subject: [PATCH 05/22] Remove Function requirement

---
 src/specializations/CylinderSurface.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl
index 9fb7d20d..3ad81737 100644
--- a/src/specializations/CylinderSurface.jl
+++ b/src/specializations/CylinderSurface.jl
@@ -9,11 +9,11 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         cyl::Meshes.CylinderSurface,
         rule::I;
         kwargs...
-) where {F <: Function, I <: IntegrationRule}
+) where {I <: IntegrationRule}
     # The generic method only parameterizes the sides
     sides = _integral(f, cyl, rule; kwargs...)
 

From 7fdf8b7c202e4b88e3b5d998607b73044af3fbad Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 11:54:46 -0500
Subject: [PATCH 06/22] Remove Function requirement

---
 src/specializations/FrustumSurface.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/specializations/FrustumSurface.jl b/src/specializations/FrustumSurface.jl
index 1de347a3..5bf1df49 100644
--- a/src/specializations/FrustumSurface.jl
+++ b/src/specializations/FrustumSurface.jl
@@ -9,11 +9,11 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         frust::Meshes.FrustumSurface,
         rule::I;
         kwargs...
-) where {F <: Function, I <: IntegrationRule}
+) where {I <: IntegrationRule}
     # The generic method only parameterizes the sides
     sides = _integral(f, frust, rule; kwargs...)
 

From 367ebb175f849ca1e30e09baa73ae46eb93b2534 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:05:41 -0500
Subject: [PATCH 07/22] Remove Function requirement

---
 src/specializations/Line.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index cb87717e..0c7da2c1 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -8,12 +8,12 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         line::Meshes.Line,
         rule::GaussLegendre;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Line, diff_method)
 
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
@@ -33,12 +33,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         line::Meshes.Line,
         rule::GaussKronrod;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Line, diff_method)
 
     # Normalize the Line s.t. line(t) is distance t from origin point
@@ -51,12 +51,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         line::Meshes.Line,
         rule::HAdaptiveCubature;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Line, diff_method)
 
     # Normalize the Line s.t. line(t) is distance t from origin point

From ae176e0a6214fd0219fb9efdd22b288ad00b44b6 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:07:19 -0500
Subject: [PATCH 08/22] Remove Function requirement

---
 src/specializations/Plane.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index b6de0e5d..ee2706cb 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -9,12 +9,12 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         plane::Meshes.Plane,
         rule::GaussLegendre;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Plane, diff_method)
 
     # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]²
@@ -40,12 +40,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         plane::Meshes.Plane,
         rule::GaussKronrod;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Plane, diff_method)
 
     # Normalize the Plane's orthogonal vectors
@@ -61,12 +61,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         plane::Meshes.Plane,
         rule::HAdaptiveCubature;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Plane, diff_method)
 
     # Normalize the Plane's orthogonal vectors
@@ -80,8 +80,8 @@ function integral(
 
     # Integrate f over the Plane
     domainunits = _units(uplane(0, 0))
-    function integrand(x::AbstractVector)
-        f(uplane(t(x[1]), t(x[2]))) * t′(x[1]) * t′(x[2]) * domainunits^2
+    function integrand(xs)
+        f(uplane(t(xs[1]), t(xs[2]))) * t′(xs[1]) * t′(xs[2]) * domainunits^2
     end
 
     # HCubature doesn't support functions that output Unitful Quantity types

From 2b521f40f956e635fa7852484b2fd3e6b69c2f8b Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:11:27 -0500
Subject: [PATCH 09/22] Remove Function requirement

---
 src/specializations/Ray.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 84d6faf2..eede812e 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -9,12 +9,12 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         ray::Meshes.Ray,
         rule::GaussLegendre;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Ray, diff_method)
 
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
@@ -24,8 +24,8 @@ function integral(
     ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v))
 
     # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞)
-    t₁(x) = FP(1 / 2) * x + FP(1 / 2)
-    t₁′(x) = FP(1 / 2)
+    t₁(x) = (1 // 2) * x + (1 // 2)
+    t₁′(x) = (1 // 2)
     t₂(x) = x / (1 - x^2)
     t₂′(x) = (1 + x^2) / (1 - x^2)^2
     t = t₂ ∘ t₁
@@ -38,12 +38,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         ray::Meshes.Ray,
         rule::GaussKronrod;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Ray, diff_method)
 
     # Normalize the Ray s.t. ray(t) is distance t from origin point
@@ -56,12 +56,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         ray::Meshes.Ray,
         rule::HAdaptiveCubature;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Ray, diff_method)
 
     # Normalize the Ray s.t. ray(t) is distance t from origin point
@@ -73,11 +73,11 @@ function integral(
 
     # Integrate f along the Ray
     domainunits = _units(ray(0))
-    integrand(x::AbstractVector) = f(ray(t(x[1]))) * t′(x[1]) * domainunits
+    integrand(xs) = f(ray(t(xs[1]))) * t′(xs[1]) * domainunits
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = zeros(FP, 1)
+    testpoint_parametriccoord = _zeros(FP, 1)
     integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))

From b1c8c39d84057b92fd3c45ef038fa2b8574db874 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:11:56 -0500
Subject: [PATCH 10/22] Remove Function requirement

---
 src/specializations/Ring.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl
index b2490ee2..36cbac64 100644
--- a/src/specializations/Ring.jl
+++ b/src/specializations/Ring.jl
@@ -27,11 +27,11 @@ calculating Jacobians that are used to calculate differential elements
 - `FP = Float64`: the floating point precision desired
 """
 function integral(
-        f::F,
+        f,
         ring::Meshes.Ring,
         rule::I;
         kwargs...
-) where {F <: Function, I <: IntegrationRule}
+) where {I <: IntegrationRule}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, rule; kwargs...), Meshes.segments(ring))
 end

From 7dcf572e00408ff4fca24d03cb904173cee263ae Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:12:17 -0500
Subject: [PATCH 11/22] Remove Function requirement

---
 src/specializations/Rope.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 7fcefa7f..34797d19 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -27,11 +27,11 @@ calculating Jacobians that are used to calculate differential elements
 - `FP = Float64`: the floating point precision desired
 """
 function integral(
-        f::F,
+        f,
         rope::Meshes.Rope,
         rule::I;
         kwargs...
-) where {F <: Function, I <: IntegrationRule}
+) where {I <: IntegrationRule}
     # Convert the Rope into Segments, sum the integrals of those
     return sum(segment -> integral(f, segment, rule; kwargs...), Meshes.segments(rope))
 end

From 560d0df82c706c73d026d3b0267070b4983d6e29 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:13:38 -0500
Subject: [PATCH 12/22] Remove Function requirement

---
 src/specializations/Tetrahedron.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index 9fb2af9a..c877397a 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -10,41 +10,41 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         tetrahedron::Meshes.Tetrahedron,
         rule::GaussLegendre;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _error_unsupported_combination("Tetrahedron", "GaussLegendre")
 end
 
 function integral(
-        f::F,
+        f,
         tetrahedron::Meshes.Tetrahedron,
         rule::GaussKronrod;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Tetrahedron, diff_method)
 
-    nil = zero(FP)
+    o = 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]
-    ∫w(w) = QuadGK.quadgk(v -> ∫vw(v, w), nil, FP(1 - w); rule.kwargs...)[1]
-    ∫ = QuadGK.quadgk(∫w, nil, one(FP); rule.kwargs...)[1]
+    ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), o, FP(1 - v - w); rule.kwargs...)[1]
+    ∫w(w) = QuadGK.quadgk(v -> ∫vw(v, w), o, FP(1 - w); rule.kwargs...)[1]
+    ∫ = QuadGK.quadgk(∫w, o, one(FP); rule.kwargs...)[1]
 
     # Apply barycentric domain correction (volume: 1/6 → actual)
     return 6 * Meshes.volume(tetrahedron) * ∫
 end
 
 function integral(
-        f::F,
+        f,
         tetrahedron::Meshes.Tetrahedron,
         rule::HAdaptiveCubature;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature")
 end
 

From 7a7e838557eee4a4207aeb764eb1edf1ddb69950 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:15:51 -0500
Subject: [PATCH 13/22] Remove Function requirement

---
 src/specializations/Triangle.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index ea0c7c11..1df8b452 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -10,12 +10,12 @@
 ################################################################################
 
 function integral(
-        f::F,
+        f,
         triangle::Meshes.Triangle,
         rule::GaussLegendre;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Triangle, diff_method)
 
     # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2
@@ -26,8 +26,8 @@ function integral(
     # Domain transformations:
     #   xᵢ [-1,1] ↦ R [0,1]
     #   xⱼ [-1,1] ↦ φ [0,π/2]
-    uR(xᵢ) = T(1 // 2) * (xᵢ + 1)
-    uφ(xⱼ) = T(π / 4) * (xⱼ + 1)
+    uR(xᵢ) = (1 // 2) * (xᵢ + 1)
+    uφ(xⱼ) = (π // 4) * (xⱼ + 1)
 
     # Integrate the Barycentric triangle by transforming it into polar coordinates
     #   with a modified radius
@@ -49,12 +49,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         triangle::Meshes.Triangle,
         rule::GaussKronrod;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Triangle, diff_method)
 
     # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0)
@@ -67,12 +67,12 @@ function integral(
 end
 
 function integral(
-        f::F,
+        f,
         triangle::Meshes.Triangle,
         rule::HAdaptiveCubature;
         diff_method::DM = Analytical(),
         FP::Type{T} = Float64
-) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat}
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _guarantee_analytical(Meshes.Triangle, diff_method)
 
     # Integrate the Barycentric triangle by transforming it into polar coordinates

From 6dc19490b74f999ccad9103e1f5b209803d07b60 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 12:46:08 -0500
Subject: [PATCH 14/22] Convert pi to FP

---
 src/specializations/Triangle.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index 1df8b452..28af5290 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -27,7 +27,7 @@ function integral(
     #   xᵢ [-1,1] ↦ R [0,1]
     #   xⱼ [-1,1] ↦ φ [0,π/2]
     uR(xᵢ) = (1 // 2) * (xᵢ + 1)
-    uφ(xⱼ) = (π // 4) * (xⱼ + 1)
+    uφ(xⱼ) = (T(π) / 4) * (xⱼ + 1)
 
     # Integrate the Barycentric triangle by transforming it into polar coordinates
     #   with a modified radius

From 5f6314bc9766712c88a0900f56b1ee9b742755ae Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 13:21:49 -0500
Subject: [PATCH 15/22] Update docstring

---
 src/integral.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/integral.jl b/src/integral.jl
index e750220f..4bb8ffc3 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -10,7 +10,7 @@ a `geometry` using a particular numerical integration `rule` with floating point
 precision of type `FP`.
 
 # Arguments
-- `f`: an integrand function with a method `f(::Meshes.Point)`
+- `f`: an integrand function, i.e. any callable 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)

From a2706d84c5ea568837e1da505d771c41df477cab Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 13:22:23 -0500
Subject: [PATCH 16/22] Update docstring

---
 src/specializations/BezierCurve.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index bd465978..f8082e1b 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -19,7 +19,7 @@
 Like [`integral`](@ref) but integrates along the domain defined by `curve`.
 
 # Arguments
-- `f`: an integrand function with a method `f(::Meshes.Point)`
+- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
 - `curve`: a `Meshes.BezierCurve` that defines the integration domain
 - `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
 

From 836032ec806985843e8e4ccff00bd3a222d6da2f Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 13:23:35 -0500
Subject: [PATCH 17/22] Update docstring [skip CI]

---
 src/specializations/Ring.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl
index 36cbac64..12d21348 100644
--- a/src/specializations/Ring.jl
+++ b/src/specializations/Ring.jl
@@ -17,7 +17,7 @@ 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)`
+- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
 - `ring`: a `Ring` that defines the integration domain
 - `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
 

From 6f7a77377f2c14a72307fb04afaa585c422459c8 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 13:24:00 -0500
Subject: [PATCH 18/22] Update docstring

---
 src/specializations/Rope.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 34797d19..e917708b 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -17,7 +17,7 @@ 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)`
+- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
 - `rope`: a `Rope` that defines the integration domain
 - `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
 

From d9108210dbb0179deee7ebda100c6124e1a3dcf1 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 14:55:32 -0500
Subject: [PATCH 19/22] Update to explain callable f

---
 README.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 1e51d48c..22ebade2 100644
--- a/README.md
+++ b/README.md
@@ -25,12 +25,12 @@ These solvers have support for integrand functions that produce scalars, vectors
 ```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.
+Performs a numerical integration of some integrand function `f` over the domain specified by `geometry`. The integrand function can be anything callable with a method `f(::Meshes.Point)`. 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()`.
+Performs a numerical integration of some integrand function `f` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. The integrand function can be anything callable with a method `f(::Meshes.Point)`.  
 
 Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics.
 

From 1904d853b404504270e31b46ad4c2b2eedf90d60 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 15:13:00 -0500
Subject: [PATCH 20/22] Update to reflect changed README

---
 docs/src/index.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index f947eb64..043d5cd1 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -25,12 +25,12 @@ These solvers have support for integrand functions that produce scalars, vectors
 ```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.
+Performs a numerical integration of some integrand function `f` over the domain specified by `geometry`. The integrand function can be anything callable with a method `f(::Meshes.Point)`. 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()`.
+Performs a numerical integration of some integrand function `f` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. The integrand function can be anything callable with a method `f(::Meshes.Point)`.
 
 Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics.
 

From 15e448d78839286cb422b8bfb77308808c9eed47 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 15:28:14 -0500
Subject: [PATCH 21/22] Add testing for callable struct as f

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

diff --git a/test/combinations.jl b/test/combinations.jl
index aa5e87f5..3d86f8e7 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -101,26 +101,29 @@ end
     a = π
     box = Box(Point(0), Point(a))
 
-    function f(p::P) where {P <: Meshes.Point}
+    struct Integrand end
+    fc = Integrand()
+
+    function (::Integrand)(p::Meshes.Point)
         t = ustrip(p.coords.x)
         sqrt(a^2 - t^2) * u"Ω/m"
     end
-    fv(p) = fill(f(p), 3)
+    fcv(p) = fill(fc(p), 3)
 
     # Scalar integrand
     sol = π * a^2 / 4 * u"Ω"
-    @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6
-    @test integral(f, box, GaussKronrod()) ≈ sol
-    @test integral(f, box, HAdaptiveCubature()) ≈ sol
+    @test integral(fc, box, GaussLegendre(100))≈sol rtol=1e-6
+    @test integral(fc, box, GaussKronrod()) ≈ sol
+    @test integral(fc, box, HAdaptiveCubature()) ≈ sol
 
     # Vector integrand
     vsol = fill(sol, 3)
-    @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6
-    @test integral(fv, box, GaussKronrod()) ≈ vsol
-    @test integral(fv, box, HAdaptiveCubature()) ≈ vsol
+    @test integral(fcv, box, GaussLegendre(100))≈vsol rtol=1e-6
+    @test integral(fcv, box, GaussKronrod()) ≈ vsol
+    @test integral(fcv, box, HAdaptiveCubature()) ≈ vsol
 
     # Integral aliases
-    @test lineintegral(f, box) ≈ sol
+    @test lineintegral(fc, box) ≈ sol
     @test_throws "not supported" surfaceintegral(f, box)
     @test_throws "not supported" volumeintegral(f, box)
 end

From 2913f807a41d1ce8df4d407cf9a1f238855c9c0c Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 28 Nov 2024 15:44:26 -0500
Subject: [PATCH 22/22] Update function name

---
 test/combinations.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/combinations.jl b/test/combinations.jl
index 3d86f8e7..4dbd16dd 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -124,8 +124,8 @@ end
 
     # Integral aliases
     @test lineintegral(fc, box) ≈ sol
-    @test_throws "not supported" surfaceintegral(f, box)
-    @test_throws "not supported" volumeintegral(f, box)
+    @test_throws "not supported" surfaceintegral(fc, box)
+    @test_throws "not supported" volumeintegral(fc, box)
 end
 
 @testitem "Meshes.Box 2D" setup=[Setup] begin