From 4c2d777960b4dc468133dc19c7a7fdd9158071bc Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 24 Nov 2024 22:53:14 -0500 Subject: [PATCH] Implement _zeros and _ones utils --- src/integral.jl | 2 +- src/specializations/BezierCurve.jl | 2 +- src/specializations/Line.jl | 4 ++-- src/specializations/Plane.jl | 2 +- src/specializations/Ray.jl | 2 +- src/specializations/Triangle.jl | 2 +- src/utils.jl | 8 ++++++++ 7 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 1a1c7878..0fee71fb 100644 --- a/src/integral.jl +++ b/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 diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index fc5905b3..7b0f44ab 100644 --- a/src/specializations/BezierCurve.jl +++ b/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 diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 4d4c4c89..88f90de1 100644 --- a/src/specializations/Line.jl +++ b/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 diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 1b34405e..b20af9da 100644 --- a/src/specializations/Plane.jl +++ b/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 diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index 9008b4a8..84d6faf2 100644 --- a/src/specializations/Ray.jl +++ b/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 diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index c63e715b..ea0c7c11 100644 --- a/src/specializations/Triangle.jl +++ b/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) .* ∫ diff --git a/src/utils.jl b/src/utils.jl index e98ade74..e69a38ef 100644 --- a/src/utils.jl +++ b/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 ################################################################################