Implement _zeros and _ones utils

src/integral.jl

@@ -100,7 +100,7 @@ function _integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP, N), ones(FP, N); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, _zeros(FP, N), _ones(FP, N); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits

src/specializations/BezierCurve.jl

@@ -82,7 +82,7 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP, 1), ones(FP, 1); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, _zeros(FP, 1), _ones(FP, 1); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits

src/specializations/Line.jl

@@ -72,12 +72,12 @@ function integral(
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = FP[0.5]
+    testpoint_parametriccoord = (FP(0.5),)
     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))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, (-one(FP),), (one(FP),); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, -_ones(FP, 1), _ones(FP, 1); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits

src/specializations/Plane.jl

@@ -91,7 +91,7 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, -ones(FP, 2), ones(FP, 2); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, -_ones(FP, 2), _ones(FP, 2); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits

src/specializations/Ray.jl

@@ -82,7 +82,7 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP, 1), ones(FP, 1); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, _zeros(FP, 1), _ones(FP, 1); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits

src/specializations/Triangle.jl

@@ -87,7 +87,7 @@ function integral(
         v = R * (1 - b / (a + b))
         return f(triangle(u, v)) * R / (a + b)^2
     end
-    ∫ = HCubature.hcubature(integrand, zeros(FP, 2), FP[1, π / 2], rule.kwargs...)[1]
+    ∫ = HCubature.hcubature(integrand, _zeros(FP, 2), (FP(1), FP(π / 2)), rule.kwargs...)[1]
 
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
     return 2 * Meshes.area(triangle) .* ∫

src/utils.jl

@@ -14,6 +14,14 @@ function _error_unsupported_combination(geometry, rule)
     throw(ArgumentError(msg))
 end
 
+# Return an NTuple{N, T} of zeros; same interface as Base.zeros() but faster
+_zeros(T::DataType, N::Int64) = ntuple(_ -> zero(T), N)
+_zeros(N::Int) = _zeros(Float64, N)
+
+# Return an NTuple{N, T} of ones; same interface as Base.ones() but faster
+_ones(T::DataType, N::Int64) = ntuple(_ -> one(T), N)
+_ones(N::Int) = _ones(Float64, N)
+
 ################################################################################
 #                           DifferentiationMethod
 ################################################################################