From 4c2d777960b4dc468133dc19c7a7fdd9158071bc Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 24 Nov 2024 22:53:14 -0500 Subject: [PATCH 1/7] 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 ################################################################################ From d4a626802897ec9cece32ae64a972516f2efff1b Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 24 Nov 2024 23:13:33 -0500 Subject: [PATCH 2/7] Change handling of negative _ones --- src/specializations/Line.jl | 2 +- src/specializations/Plane.jl | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 88f90de1..4bebeb1d 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -77,7 +77,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, 1), _ones(FP, 1); rule.kwargs...)[1] + value = HCubature.hcubature(uintegrand, (-_one(FP),), (_one(FP),); rule.kwargs...)[1] # Reapply units return value .* integrandunits diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index b20af9da..685ca7b4 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -91,7 +91,9 @@ 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] + a = 0 .- _ones(FP, 2) + b = _ones(FP, 2) + value = HCubature.hcubature(uintegrand, a, b; rule.kwargs...)[1] # Reapply units return value .* integrandunits From 921b18627c109fcc56b918fa8c0ac9a158a19cf8 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 24 Nov 2024 23:35:16 -0500 Subject: [PATCH 3/7] Generalize integrand arg --- src/specializations/Line.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 4bebeb1d..7e7a3d91 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -68,7 +68,7 @@ function integral( # Integrate f along the Line differential(line, x) = t′(x) * _units(line(0)) - integrand(x::AbstractVector) = f(line(t(x[1]))) * differential(line, x[1]) + integrand(xs) = f(line(t(xs[1]))) * differential(line, xs[1]) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f From 8f9483ab292aa0a4500b74e39e343662b3139eb0 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 24 Nov 2024 23:51:44 -0500 Subject: [PATCH 4/7] Fix typo --- src/specializations/Line.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 7e7a3d91..cb87717e 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -77,7 +77,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, (-_one(FP),), (_one(FP),); rule.kwargs...)[1] + value = HCubature.hcubature(uintegrand, (-one(FP),), (one(FP),); rule.kwargs...)[1] # Reapply units return value .* integrandunits From 7dde5721c4d2b3624bdb95f1ca0e8e67ac1be531 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 25 Nov 2024 00:06:38 -0500 Subject: [PATCH 5/7] Add new tests --- test/utils.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/utils.jl b/test/utils.jl index f5c26308..a42ae24f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,5 +1,6 @@ @testitem "Utilities" setup=[Setup] begin using LinearAlgebra: norm + using MeshIntegrals: _units, _zeros, _ones # _KVector v = Meshes.Vec(3, 4) @@ -7,7 +8,15 @@ # _units p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") - @test MeshIntegrals._units(p) == u"m" + @test _units(p) == u"m" + + # _zeros + @test _zeros(2) == (0.0, 0.0) + @test _zeros(Float32, 2) == (0.0f0, 0.0f0) + + # _ones + @test _ones(2) == (1.0, 1.0) + @test _ones(Float32, 2) == (1.0f0, 1.0f0) end @testitem "DifferentiationMethod" setup=[Setup] begin From a57c2549dbba13de1ce1463fe5b50c1e85f21892 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 25 Nov 2024 04:45:09 -0500 Subject: [PATCH 6/7] Drop leaving zero Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/specializations/Plane.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 685ca7b4..53a2666b 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 - a = 0 .- _ones(FP, 2) + a = .- _ones(FP, 2) b = _ones(FP, 2) value = HCubature.hcubature(uintegrand, a, b; rule.kwargs...)[1] From f1952e3019cef782c9c12755b37262fc04abc2a5 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 25 Nov 2024 04:57:06 -0500 Subject: [PATCH 7/7] Apply bot formatting suggestion Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/specializations/Plane.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 53a2666b..b6de0e5d 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 - a = .- _ones(FP, 2) + a = .-_ones(FP, 2) b = _ones(FP, 2) value = HCubature.hcubature(uintegrand, a, b; rule.kwargs...)[1]