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.
 
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.
 
diff --git a/src/integral.jl b/src/integral.jl
index de69db5d..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)
@@ -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
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...
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index c3540cbc..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
 
@@ -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
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...)
 
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...)
 
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...)
 
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
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
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))
diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl
index b2490ee2..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
 
@@ -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
index 7fcefa7f..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
 
@@ -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
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
 
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index ea0c7c11..28af5290 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ⱼ) = (T(π) / 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
diff --git a/test/combinations.jl b/test/combinations.jl
index aa5e87f5..4dbd16dd 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -101,28 +101,31 @@ 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_throws "not supported" surfaceintegral(f, box)
-    @test_throws "not supported" volumeintegral(f, box)
+    @test lineintegral(fc, box) ≈ sol
+    @test_throws "not supported" surfaceintegral(fc, box)
+    @test_throws "not supported" volumeintegral(fc, box)
 end
 
 @testitem "Meshes.Box 2D" setup=[Setup] begin