Generalize integrand f to any callable (#137)

mikeingold · web-flow · commit 3635a37940a2 · 2024-11-29T10:03:53.000-05:00
* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Remove Function requirement

* Convert pi to FP

* Update docstring

* Update docstring

* Update docstring [skip CI]

* Update docstring

* Update to explain callable f

* Update to reflect changed README

* Add testing for callable struct as f

* Update function name
diff --git 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.
 
diff --git 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.
 
diff --git 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)
@@ -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...
diff --git 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...
diff --git 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
 
@@ -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)
 
diff --git 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...)
 
diff --git 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...)
 
diff --git 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...)
 
diff --git 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
diff --git 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
diff --git 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))
diff --git 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
 
@@ -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
diff --git 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
 
@@ -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
diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
diff --git a/test/combinations.jl b/test/combinations.jl
