Generalize integrand f to any callable

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.
 

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.
 

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...

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...

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)
 

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...)
 

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...)
 

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...)
 

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

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

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))

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

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