From 1f1a10fc7d2e1028fbcb31db126ef1d5b042a78b Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 15:08:58 -0500 Subject: [PATCH 01/47] Initial attempt at new test generation --- test/combinations.jl | 135 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 117 insertions(+), 18 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 2ffb1b9a..14b72c1f 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -2,33 +2,132 @@ # - All supported combinations of integral(f, ::Geometry, ::IntegrationAlgorithm) produce accurate results # - Invalid applications of integral aliases (e.g. lineintegral) produce a descriptive error -@testitem "Meshes.Ball 2D" setup=[Setup] begin +#=============================================================================== + Test Generation Infrastructure +===============================================================================# + +@testsnippet Combinations begin + using LinearAlgebra: norm + using Meshes + using MeshIntegrals + using Unitful + + struct Callable{F <: Function} + f::F + end + (c::Callable)(p) = c.f(p) + + # Stores a testable combination + struct TestableGeometry{F <: Function, G <: Geometry, U <: Unitful.Quantity} + integrand::F + geometry::G + solution::U + end + + # Indicates which functions/rules are supported for a particular geometry + struct SupportStatus + lineintegral::Bool + surfaceintegral::Bool + volumeintegral::Bool + gausskronrod::Bool + gausslegendre::Bool + hadaptivecubature::Bool + end + + # Shortcut constructor for geometries with typical support structure + function SupportStatus(sym::Symbol) + if sym == :line + aliases = Bool.((0, 0, 1)) + rules = Bool.((1, 1, 1)) + return SupportStatus(aliases..., rules...) + elseif sym == :surface + aliases = Bool.((0, 1, 0)) + rules = Bool.((1, 1, 1)) + return SupportStatus(aliases..., rules...) + elseif sym == :volume + aliases = Bool.((0, 0, 1)) + rules = Bool.((0, 1, 1)) + return SupportStatus(aliases..., rules...) + else + error("Unrecognized SupportStatus shortcut $(string(sym))") + end + end + + function runtests(testable::TestableGeometry, supports::SupportStatus) + #= + for fname in ("lineintegral", "surfaceintegral", "volumeintegral") + if supports.$fname + @test $fname(testable.integrand, testable.geometry) ≈ testable.solution + else + @test_throws "not supported" $fname(testable.integrand, testable.geometry) + end + end + =# + + # lineintegral + if supports.lineintegral + @test lineintegral(testable.integrand, testable.geometry) ≈ testable.solution + else + @test_throws "not supported" lineintegral(testable.integrand, testable.geometry) + end + + # surfaceintegral + if supports.surfaceintegral + @test surfaceintegral(testable.integrand, testable.geometry) ≈ testable.solution + else + @test_throws "not supported" surfaceintegral(testable.integrand, testable.geometry) + end + + # volumeintegral + if supports.volumeintegral + @test volumeintegral(testable.integrand, testable.geometry) ≈ testable.solution + else + @test_throws "not supported" volumeintegral(testable.integrand, testable.geometry) + end + + if supports.gausskronrod + # Scalar integrand + sol = testable.solution + @test integral(testable.integrand, testable.geometry, GaussKronrod()) ≈ sol + + # Vector integrand + fv(p) = fill(testable.integrand(p), 3) + sol_v = fill(testable.solution, 3) + @test integral(fv, testable.geometry, GaussKronrod()) ≈ sol_v + + # Callable integrand + f = Callable(testable.integrand) + @test integral(f, testable.geometry, GaussKronrod()) ≈ sol + else + @test_throws "not supported" integral(testable.integrand, testable.geometry, GaussKronrod()) + end + + # supports.gausslegendre, GaussLegendre(100) + + # supports.hadaptivecubature, HAdaptiveCubature() + end + +end + +#=============================================================================== + Create and Test Geometries +===============================================================================# + +@testitem "Meshes.Ball 2D" setup=[Combinations] begin origin = Point(0, 0) radius = 2.8 ball = Ball(origin, radius) - function f(p::P) where {P <: Meshes.Point} + function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) exp(-r^2) end - fv(p) = fill(f(p), 3) - # Scalar integrand - sol = (π - π * exp(-radius^2)) * u"m^2" - @test integral(f, ball, GaussLegendre(100)) ≈ sol - @test integral(f, ball, GaussKronrod()) ≈ sol - @test integral(f, ball, HAdaptiveCubature()) ≈ sol + solution = (π - π * exp(-radius^2)) * u"m^2" - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, ball, GaussLegendre(100)) ≈ vsol - @test integral(fv, ball, GaussKronrod()) ≈ vsol - @test integral(fv, ball, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, ball) - @test surfaceintegral(f, ball) ≈ sol - @test_throws "not supported" volumeintegral(f, ball) + support = SupportStatus(:surface) + testable = TestableGeometry(integrand, ball, solution) + runtests(testable, support) end @testitem "Meshes.Ball 3D" setup=[Setup] begin From 339510208d4bc2f8f3356ad7e2ddd654440ab792 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 17:01:09 -0500 Subject: [PATCH 02/47] Generate alias tests in a loop --- test/combinations.jl | 32 +++++--------------------------- 1 file changed, 5 insertions(+), 27 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 14b72c1f..f881472a 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -54,36 +54,14 @@ end function runtests(testable::TestableGeometry, supports::SupportStatus) - #= - for fname in ("lineintegral", "surfaceintegral", "volumeintegral") - if supports.$fname - @test $fname(testable.integrand, testable.geometry) ≈ testable.solution + for alias in (lineintegral, surfaceintegral, volumeintegral) + alias_symbol = first(methods(alias)).name + if getfield(supports, alias_symbol) + @test alias(testable.integrand, testable.geometry) ≈ testable.solution else - @test_throws "not supported" $fname(testable.integrand, testable.geometry) + @test_throws "not supported" alias(testable.integrand, testable.geometry) end end - =# - - # lineintegral - if supports.lineintegral - @test lineintegral(testable.integrand, testable.geometry) ≈ testable.solution - else - @test_throws "not supported" lineintegral(testable.integrand, testable.geometry) - end - - # surfaceintegral - if supports.surfaceintegral - @test surfaceintegral(testable.integrand, testable.geometry) ≈ testable.solution - else - @test_throws "not supported" surfaceintegral(testable.integrand, testable.geometry) - end - - # volumeintegral - if supports.volumeintegral - @test volumeintegral(testable.integrand, testable.geometry) ≈ testable.solution - else - @test_throws "not supported" volumeintegral(testable.integrand, testable.geometry) - end if supports.gausskronrod # Scalar integrand From a4c5a5d970b200ccf423567dd1f66cabf15211da Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:16:51 -0500 Subject: [PATCH 03/47] Loop over all integration rules --- test/combinations.jl | 68 ++++++++++++++++++++------------------------ 1 file changed, 31 insertions(+), 37 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index f881472a..8bb07de3 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -63,26 +63,31 @@ end end - if supports.gausskronrod - # Scalar integrand - sol = testable.solution - @test integral(testable.integrand, testable.geometry, GaussKronrod()) ≈ sol - - # Vector integrand - fv(p) = fill(testable.integrand(p), 3) - sol_v = fill(testable.solution, 3) - @test integral(fv, testable.geometry, GaussKronrod()) ≈ sol_v - - # Callable integrand - f = Callable(testable.integrand) - @test integral(f, testable.geometry, GaussKronrod()) ≈ sol - else - @test_throws "not supported" integral(testable.integrand, testable.geometry, GaussKronrod()) + iter_rules = ( + (supports.gausskronrod, GaussKronrod()), + (supports.gausslegendre, GaussLegendre(100)), + (supports.hadaptivecubature, HAdaptiveCubature()) + ) + + for (supported, rule) in iter_rules + if supported + # Scalar integrand + sol = testable.solution + @test integral(testable.integrand, testable.geometry, rule) ≈ sol + + # Callable integrand + f = Callable(testable.integrand) + @test integral(f, testable.geometry, rule) ≈ sol + + # Vector integrand + fv(p) = fill(testable.integrand(p), 3) + sol_v = fill(testable.solution, 3) + @test integral(fv, testable.geometry, rule) ≈ sol_v + else + @test_throws "not supported" integral(testable.integrand, testable.geometry, rule) + end end - # supports.gausslegendre, GaussLegendre(100) - - # supports.hadaptivecubature, HAdaptiveCubature() end end @@ -108,38 +113,27 @@ end runtests(testable, support) end -@testitem "Meshes.Ball 3D" setup=[Setup] begin +@testitem "Meshes.Ball 3D" setup=[Combinations] begin using SpecialFunctions: erf center = Point(1, 2, 3) radius = 2.8u"m" ball = Ball(center, radius) - function f(p::P) where {P <: Meshes.Point} + function integrand(p::P) where {P <: Meshes.Point} offset = p - center ur = norm(offset) r = ustrip(u"m", ur) exp(-r^2) end - fv(p) = fill(f(p), 3) - # Scalar integrand - r = ustrip(u"m", radius) - sol = (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" - @test integral(f, ball, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, ball, GaussKronrod())≈sol - @test integral(f, ball, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, ball, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, ball, GaussKronrod())≈vsol - @test integral(fv, ball, HAdaptiveCubature()) ≈ vsol + solution = let r = ustrip(u"m", radius) + (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" + end - # Integral aliases - @test_throws "not supported" lineintegral(f, ball) - @test_throws "not supported" surfaceintegral(f, ball) - @test volumeintegral(f, ball) ≈ sol + support = SupportStatus(:volume) + testable = TestableGeometry(integrand, ball, solution) + runtests(testable, support) end @testitem "Meshes.BezierCurve" setup=[Setup] begin From 9620a2d5a808a115cc2b300dfa7099a5ff204c2d Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:35:57 -0500 Subject: [PATCH 04/47] Implement for Circle --- test/combinations.jl | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 8bb07de3..a445f985 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -285,36 +285,27 @@ end @test_throws "not supported" volumeintegral(f, box) end -@testitem "Meshes.Circle" setup=[Setup] begin +@testitem "Meshes.Circle" setup=[Combinations] begin + # Geometry center = Point(1, 2, 3) n̂ = Vec(1 / 2, 1 / 2, sqrt(2) / 2) plane = Plane(center, n̂) radius = 4.4 circle = Circle(plane, radius) - function f(p::P) where {P <: Meshes.Point} + # Integrand + function integrand(p::P) where {P <: Meshes.Point} offset = p - center r = ustrip(u"m", norm(offset)) exp(-r^2) end - fv(p) = fill(f(p), 3) # Scalar integrand - sol = 2π * radius * exp(-radius^2) * u"m" - @test integral(f, circle, GaussLegendre(100)) ≈ sol - @test integral(f, circle, GaussKronrod()) ≈ sol - @test integral(f, circle, HAdaptiveCubature()) ≈ sol + solution = 2π * radius * exp(-radius^2) * u"m" - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, circle, GaussLegendre(100)) ≈ vsol - @test integral(fv, circle, GaussKronrod()) ≈ vsol - @test integral(fv, circle, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test lineintegral(f, circle) ≈ sol - @test_throws "not supported" surfaceintegral(f, circle) - @test_throws "not supported" volumeintegral(f, circle) + # Package and run tests + testable = TestableGeometry(integrand, circle, solution) + runtests(testable, SupportStatus(:line)) end @testitem "Meshes.Cone" setup=[Setup] begin From 49e983b1e0206c146a37bdec8d8ab2f3c6783171 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:36:57 -0500 Subject: [PATCH 05/47] Code cleanup --- test/combinations.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index a445f985..aae57d08 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -1,6 +1,6 @@ -# This section tests for: -# - All supported combinations of integral(f, ::Geometry, ::IntegrationAlgorithm) produce accurate results -# - Invalid applications of integral aliases (e.g. lineintegral) produce a descriptive error +# This section tests: +# - All supported combinations of integral(f, ::Geometry, ::IntegrationAlgorithm) +# - Invalid applications of integral aliases produce a descriptive error #=============================================================================== Test Generation Infrastructure @@ -54,6 +54,7 @@ end function runtests(testable::TestableGeometry, supports::SupportStatus) + # Test alias functions for alias in (lineintegral, surfaceintegral, volumeintegral) alias_symbol = first(methods(alias)).name if getfield(supports, alias_symbol) @@ -69,6 +70,7 @@ (supports.hadaptivecubature, HAdaptiveCubature()) ) + # Test rules for (supported, rule) in iter_rules if supported # Scalar integrand @@ -97,29 +99,34 @@ end ===============================================================================# @testitem "Meshes.Ball 2D" setup=[Combinations] begin + # Geometry origin = Point(0, 0) radius = 2.8 ball = Ball(origin, radius) + # Integrand function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) exp(-r^2) end + # Solution solution = (π - π * exp(-radius^2)) * u"m^2" - support = SupportStatus(:surface) + # Package and run tests testable = TestableGeometry(integrand, ball, solution) - runtests(testable, support) + runtests(testable, SupportStatus(:surface)) end @testitem "Meshes.Ball 3D" setup=[Combinations] begin using SpecialFunctions: erf + # Geometry center = Point(1, 2, 3) radius = 2.8u"m" ball = Ball(center, radius) + # Integrand function integrand(p::P) where {P <: Meshes.Point} offset = p - center ur = norm(offset) @@ -127,18 +134,21 @@ end exp(-r^2) end + # Solution solution = let r = ustrip(u"m", radius) (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" end - support = SupportStatus(:volume) + # Package and run tests testable = TestableGeometry(integrand, ball, solution) - runtests(testable, support) + runtests(testable, SupportStatus(:volume)) end @testitem "Meshes.BezierCurve" setup=[Setup] begin - curve = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)]) + # Geometry + curve = BezierCurve([Point(t, sin(t), 0) for t in range(-π, π, length = 361)]) + # Integrand function f(p::P) where {P <: Meshes.Point} ux = ustrip(p.coords.x) (1 / sqrt(1 + cos(ux)^2)) * u"Ω/m" From d9b20dfbbeabfe770727707e8750939d80e68c68 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:37:25 -0500 Subject: [PATCH 06/47] Bugfix --- test/combinations.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index aae57d08..5ad8399b 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -37,7 +37,7 @@ # Shortcut constructor for geometries with typical support structure function SupportStatus(sym::Symbol) if sym == :line - aliases = Bool.((0, 0, 1)) + aliases = Bool.((1, 0, 0)) rules = Bool.((1, 1, 1)) return SupportStatus(aliases..., rules...) elseif sym == :surface @@ -89,7 +89,6 @@ @test_throws "not supported" integral(testable.integrand, testable.geometry, rule) end end - end end From b9c0fd4ea0b8af4e16eea825806a6885ed351d45 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:41:29 -0500 Subject: [PATCH 07/47] Implement for Cone --- test/combinations.jl | 31 ++++++++++--------------------- 1 file changed, 10 insertions(+), 21 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 5ad8399b..c2a5cdd3 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -317,36 +317,25 @@ end runtests(testable, SupportStatus(:line)) end -@testitem "Meshes.Cone" setup=[Setup] begin +@testitem "Meshes.Cone" setup=[Combinations] begin + # Geometry r = 2.5u"m" - h = 2.5u"m" + h = 3.5u"m" origin = Point(0, 0, 0) xy_plane = Plane(origin, Vec(0, 0, 1)) base = Disk(xy_plane, r) apex = Point(0.0u"m", 0.0u"m", h) cone = Cone(base, apex) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - - _volume_cone_rightcircular(h, r) = π * r^2 * h / 3 - - # Scalar integrand - sol = _volume_cone_rightcircular(r, h) - @test integral(f, cone, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, cone, GaussKronrod()) - @test integral(f, cone, HAdaptiveCubature()) ≈ sol + # Integrand + integrand(p) = 1.0u"A" - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, cone, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, cone, GaussKronrod()) - @test integral(fv, cone, HAdaptiveCubature()) ≈ vsol + # Solution + solution = (π * r^2 * h / 3) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, cone) - @test_throws "not supported" surfaceintegral(f, cone) - @test volumeintegral(f, cone) ≈ sol + # Package and run tests + testable = TestableGeometry(integrand, cone, solution) + runtests(testable, SupportStatus(:volume)) end @testitem "Meshes.ConeSurface" setup=[Setup] begin From 6bccee1c7564111e42eccf0e068c1a602189a2c8 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:45:23 -0500 Subject: [PATCH 08/47] Implement for Cylinder --- test/combinations.jl | 27 +++++++++------------------ 1 file changed, 9 insertions(+), 18 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index c2a5cdd3..294239ec 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -370,30 +370,21 @@ end @test_throws "not supported" volumeintegral(f, cone) end -@testitem "Meshes.Cylinder" setup=[Setup] begin +@testitem "Meshes.Cylinder" setup=[Combinations] begin + # Geometry pt_w = Point(-1, 0, 0) pt_e = Point(1, 0, 0) cyl = Cylinder(pt_e, pt_w, 2.5) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = Meshes.measure(cyl) - @test integral(f, cyl, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, cyl, GaussKronrod()) - @test integral(f, cyl, HAdaptiveCubature()) ≈ sol + # Integrand + integrand(p) = 1.0u"A" - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, cyl, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, cyl, GaussKronrod()) - @test integral(fv, cyl, HAdaptiveCubature()) ≈ vsol + # Solution + solution = Meshes.measure(cyl) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, cyl) - @test_throws "not supported" surfaceintegral(f, cyl) - @test volumeintegral(f, cyl) ≈ sol + # Package and run tests + testable = TestableGeometry(integrand, cyl, solution) + runtests(testable, SupportStatus(:volume)) end @testitem "Meshes.CylinderSurface" setup=[Setup] begin From fe5138ced17660bbb6eece3c2003ea0abf851fd9 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:52:10 -0500 Subject: [PATCH 09/47] Implement for CylinderSurface and Disk --- test/combinations.jl | 60 +++++++++++++------------------------------- 1 file changed, 18 insertions(+), 42 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 294239ec..0de6148c 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -376,10 +376,8 @@ end pt_e = Point(1, 0, 0) cyl = Cylinder(pt_e, pt_w, 2.5) - # Integrand + # Integrand & Solution integrand(p) = 1.0u"A" - - # Solution solution = Meshes.measure(cyl) * u"A" # Package and run tests @@ -387,62 +385,40 @@ end runtests(testable, SupportStatus(:volume)) end -@testitem "Meshes.CylinderSurface" setup=[Setup] begin +@testitem "Meshes.CylinderSurface" setup=[Combinations] begin + # Geometry pt_w = Point(-1, 0, 0) pt_e = Point(1, 0, 0) cyl = CylinderSurface(pt_e, pt_w, 2.5) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = Meshes.measure(cyl) - @test integral(f, cyl, GaussLegendre(100)) ≈ sol - @test integral(f, cyl, GaussKronrod()) ≈ sol - @test integral(f, cyl, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, cyl, GaussLegendre(100)) ≈ vsol - @test integral(fv, cyl, GaussKronrod()) ≈ vsol - @test integral(fv, cyl, HAdaptiveCubature()) ≈ vsol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(cyl) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, cyl) - @test surfaceintegral(f, cyl) ≈ sol - @test_throws "not supported" volumeintegral(f, cyl) + # Package and run tests + testable = TestableGeometry(integrand, cyl, solution) + runtests(testable, SupportStatus(:surface)) end -@testitem "Meshes.Disk" setup=[Setup] begin +@testitem "Meshes.Disk" setup=[Combinations] begin + # Geometry center = Point(1, 2, 3) n̂ = Vec(1 / 2, 1 / 2, sqrt(2) / 2) plane = Plane(center, n̂) radius = 2.5 disk = Disk(plane, radius) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} offset = p - center r = ustrip(u"m", norm(offset)) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) + solution = (π - π * exp(-radius^2)) * u"A*m^2" - # Scalar integrand - sol = (π - π * exp(-radius^2)) * u"m^2" - @test integral(f, disk, GaussLegendre(100)) ≈ sol - @test integral(f, disk, GaussKronrod()) ≈ sol - @test integral(f, disk, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, disk, GaussLegendre(100)) ≈ vsol - @test integral(fv, disk, GaussKronrod()) ≈ vsol - @test integral(fv, disk, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, disk) - @test surfaceintegral(f, disk) ≈ sol - @test_throws "not supported" volumeintegral(f, disk) + # Package and run tests + testable = TestableGeometry(integrand, disk, solution) + runtests(testable, SupportStatus(:surface)) end @testitem "Meshes.Ellipsoid" setup=[Setup] begin From 9ab76007df316ec08c9276eb88f26625374fb19f Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 21:56:17 -0500 Subject: [PATCH 10/47] Formatting --- test/combinations.jl | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 0de6148c..9a9b1ad0 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -39,7 +39,7 @@ if sym == :line aliases = Bool.((1, 0, 0)) rules = Bool.((1, 1, 1)) - return SupportStatus(aliases..., rules...) + return SupportStatus(aliases..., rules...) elseif sym == :surface aliases = Bool.((0, 1, 0)) rules = Bool.((1, 1, 1)) @@ -87,11 +87,10 @@ @test integral(fv, testable.geometry, rule) ≈ sol_v else @test_throws "not supported" integral(testable.integrand, testable.geometry, rule) - end - end - end - -end + end # if + end # for + end # function +end #testsnippet #=============================================================================== Create and Test Geometries @@ -103,13 +102,11 @@ end radius = 2.8 ball = Ball(origin, radius) - # Integrand + # Integrand & Solution function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) exp(-r^2) end - - # Solution solution = (π - π * exp(-radius^2)) * u"m^2" # Package and run tests @@ -125,15 +122,13 @@ end radius = 2.8u"m" ball = Ball(center, radius) - # Integrand + # Integrand & Solution function integrand(p::P) where {P <: Meshes.Point} offset = p - center ur = norm(offset) r = ustrip(u"m", ur) exp(-r^2) end - - # Solution solution = let r = ustrip(u"m", radius) (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" end From b2147850dbb5bd666bd2799f63b17a39d8ee2c32 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 22:04:26 -0500 Subject: [PATCH 11/47] Implement for Line, ParaboloidSurface, and ParametrizedCurve --- test/combinations.jl | 97 +++++++++++++------------------------------- 1 file changed, 28 insertions(+), 69 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 9a9b1ad0..13fd9c53 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -508,103 +508,62 @@ end end @testitem "Meshes.Line" setup=[Setup] begin - a = Point(0.0u"m", 0.0u"m", 0.0u"m") - b = Point(1.0u"m", 1.0u"m", 1.0u"m") + # Geometry + a = Point(0, 0, 0) + b = Point(1, 1, 1) line = Line(a, b) + # Integrand & solution function f(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = sqrt(π) * u"m" - @test integral(f, line, GaussLegendre(100)) ≈ sol - @test integral(f, line, GaussKronrod()) ≈ sol - @test integral(f, line, HAdaptiveCubature()) ≈ sol + solution = sqrt(π) * u"A*m" - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, line, GaussLegendre(100)) ≈ vsol - @test integral(fv, line, GaussKronrod()) ≈ vsol - @test integral(fv, line, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test lineintegral(f, line) ≈ sol - @test_throws "not supported" surfaceintegral(f, line) - @test_throws "not supported" volumeintegral(f, line) + # Package and run tests + testable = TestableGeometry(integrand, line, solution) + runtests(testable, SupportStatus(:line)) end -@testitem "Meshes.ParaboloidSurface" setup=[Setup] begin +@testitem "Meshes.ParaboloidSurface" setup=[Combinations] begin origin = Point(0, 0, 0) parab = ParaboloidSurface(origin, 2.5, 4.15) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = Meshes.measure(parab) - @test integral(f, parab, GaussLegendre(100)) ≈ sol - @test integral(f, parab, GaussKronrod()) ≈ sol - @test integral(f, parab, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, parab, GaussLegendre(100)) ≈ vsol - @test integral(fv, parab, GaussKronrod()) ≈ vsol - @test integral(fv, parab, HAdaptiveCubature()) ≈ vsol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(parab) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, parab) - @test surfaceintegral(f, parab) ≈ sol - @test_throws "not supported" volumeintegral(f, parab) + # Package and run tests + testable = TestableGeometry(integrand, parab, solution) + runtests(testable, SupportStatus(:surface)) end -@testitem "ParametrizedCurve" setup=[Setup] begin +@testitem "ParametrizedCurve" setup=[Combinations] begin # ParametrizedCurve has been added in Meshes v0.51.20 # If the version is specified as minimal compat bound in the Project.toml, the downgrade test fails if pkgversion(Meshes) >= v"0.51.20" using CoordRefSystems: Polar + # Geometries # Parameterize a circle centered on origin with specified radius radius = 4.4 curve_cart = ParametrizedCurve( t -> Point(radius * cos(t), radius * sin(t)), (0.0, 2π)) curve_polar = ParametrizedCurve(t -> Point(Polar(radius, t)), (0.0, 2π)) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} ur = norm(to(p)) r = ustrip(u"m", ur) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = 2π * radius * exp(-radius^2) * u"m" - @test integral(f, curve_cart, GaussLegendre(100)) ≈ sol - @test integral(f, curve_cart, GaussKronrod()) ≈ sol - @test integral(f, curve_cart, HAdaptiveCubature()) ≈ sol - @test integral(f, curve_polar, GaussLegendre(100)) ≈ sol - @test integral(f, curve_polar, GaussKronrod()) ≈ sol - @test integral(f, curve_polar, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, curve_cart, GaussLegendre(100)) ≈ vsol - @test integral(fv, curve_cart, GaussKronrod()) ≈ vsol - @test integral(fv, curve_cart, HAdaptiveCubature()) ≈ vsol - @test integral(fv, curve_polar, GaussLegendre(100)) ≈ vsol - @test integral(fv, curve_polar, GaussKronrod()) ≈ vsol - @test integral(fv, curve_polar, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test lineintegral(f, curve_cart) ≈ sol - @test_throws "not supported" surfaceintegral(f, curve_cart) - @test_throws "not supported" volumeintegral(f, curve_cart) - @test lineintegral(f, curve_polar) ≈ sol - @test_throws "not supported" surfaceintegral(f, curve_polar) - @test_throws "not supported" volumeintegral(f, curve_polar) + solution = 2π * radius * exp(-radius^2) * u"A*m" + + # Package and run tests + testable_cart = TestableGeometry(integrand, curve_cart, solution) + runtests(testable_cart, SupportStatus(:line)) + testable_polar = TestableGeometry(integrand, curve_cart, solution) + runtests(testable_polar, SupportStatus(:line)) end end From dd7d15616f3982c065eaf353d44d355a988a977a Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 4 Dec 2024 22:05:58 -0500 Subject: [PATCH 12/47] Update spelling --- test/combinations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index 13fd9c53..f124b2bd 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -545,7 +545,7 @@ end using CoordRefSystems: Polar # Geometries - # Parameterize a circle centered on origin with specified radius + # Parametrize a circle centered on origin with specified radius radius = 4.4 curve_cart = ParametrizedCurve( t -> Point(radius * cos(t), radius * sin(t)), (0.0, 2π)) From 26b1e84f1bb0b2978cfa6c268783090f5ced359c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 10:37:57 -0500 Subject: [PATCH 13/47] Update setup --- test/combinations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index f124b2bd..f8e7142f 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -507,7 +507,7 @@ end @test volumeintegral(f, hexahedron) ≈ sol end -@testitem "Meshes.Line" setup=[Setup] begin +@testitem "Meshes.Line" setup=[Combinations] begin # Geometry a = Point(0, 0, 0) b = Point(1, 1, 1) From ebb3b9be0e2c8d321067f2956169e014ec488a89 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 11:33:16 -0500 Subject: [PATCH 14/47] Implement for Plane, Quadrangle, and Ray --- test/combinations.jl | 92 +++++++++++++++----------------------------- 1 file changed, 30 insertions(+), 62 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index f8e7142f..0144211d 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -567,90 +567,58 @@ end end end -@testitem "Meshes.Plane" setup=[Setup] begin +@testitem "Meshes.Plane" setup=[Combinations] begin + # Geometry p = Point(0.0u"m", 0.0u"m", 0.0u"m") v = Vec(0.0u"m", 0.0u"m", 1.0u"m") plane = Plane(p, v) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = π * u"m^2" - @test integral(f, plane, GaussLegendre(100)) ≈ sol - @test integral(f, plane, GaussKronrod()) ≈ sol - @test integral(f, plane, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, plane, GaussLegendre(100)) ≈ vsol - @test integral(fv, plane, GaussKronrod()) ≈ vsol - @test integral(fv, plane, HAdaptiveCubature()) ≈ vsol + solution = π * u"A*m^2" - # Integral aliases - @test_throws "not supported" lineintegral(f, plane) - @test surfaceintegral(f, plane) ≈ sol - @test_throws "not supported" volumeintegral(f, plane) + # Package and run tests + testable = TestableGeometry(integrand, plane, solution) + runtests(testable, SupportStatus(:surface)) end -@testitem "Meshes.Quadrangle" setup=[Setup] begin +@testitem "Meshes.Quadrangle" setup=[Combinations] begin using SpecialFunctions: erf + + # Geometry quadrangle = Quadrangle((-1.0, 0.0), (-1.0, 1.0), (1.0, 1.0), (1.0, 0.0)) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) + solution = 0.5 * π * erf(1)^2 * u"A*m^2" - # Scalar integrand - sol = 0.5 * π * erf(1)^2 * u"m^2" - @test integral(f, quadrangle, GaussLegendre(100)) ≈ sol - @test integral(f, quadrangle, GaussKronrod()) ≈ sol - @test integral(f, quadrangle, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, quadrangle, GaussLegendre(100)) ≈ vsol - @test integral(fv, quadrangle, GaussKronrod()) ≈ vsol - @test integral(fv, quadrangle, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, quadrangle) - @test surfaceintegral(f, quadrangle) ≈ sol - @test_throws "not supported" volumeintegral(f, quadrangle) + # Package and run tests + testable = TestableGeometry(integrand, quadrangle, solution) + runtests(testable, SupportStatus(:surface)) end -@testitem "Meshes.Ray" setup=[Setup] begin - a = Point(0.0u"m", 0.0u"m", 0.0u"m") - v = Vec(1.0u"m", 1.0u"m", 1.0u"m") +@testitem "Meshes.Ray" setup=[Combinations] begin + # Geometry + a = Point(0, 0, 0) + v = Vec(1, 1, 1) ray = Ray(a, v) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) + solution = sqrt(π) / 2 * u"A*m" - # Scalar integrand - sol = sqrt(π) / 2 * u"m" - @test integral(f, ray, GaussLegendre(100)) ≈ sol - @test integral(f, ray, GaussKronrod()) ≈ sol - @test integral(f, ray, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, ray, GaussLegendre(100)) ≈ vsol - @test integral(fv, ray, GaussKronrod()) ≈ vsol - @test integral(fv, ray, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test lineintegral(f, ray) ≈ sol - @test_throws "not supported" surfaceintegral(f, ray) - @test_throws "not supported" volumeintegral(f, ray) + # Package and run tests + testable = TestableGeometry(integrand, ray, solution) + runtests(testable, SupportStatus(:line)) end @testitem "Meshes.Ring" setup=[Setup] begin From 76ef87b327f307e7aa5fd88523e0e2a2e17a9778 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 11:34:21 -0500 Subject: [PATCH 15/47] Update integrand name --- test/combinations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index 0144211d..94542b8e 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -514,7 +514,7 @@ end line = Line(a, b) # Integrand & solution - function f(p::P) where {P <: Meshes.Point} + function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end From 77de4085629701aa4cffd33c39ffa24315b96663 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 11:47:25 -0500 Subject: [PATCH 16/47] Implement for Triangle --- test/combinations.jl | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 94542b8e..bd2d5f6b 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -833,29 +833,18 @@ end @test_throws "not supported" volumeintegral(f, torus) end -@testitem "Meshes.Triangle" setup=[Setup] begin +@testitem "Meshes.Triangle" setup=[Combinations] begin + # Geometry pt_n = Point(0, 1, 0) pt_w = Point(-1, 0, 0) pt_e = Point(1, 0, 0) triangle = Triangle(pt_e, pt_n, pt_w) - f(p) = 1.0u"A" - fv(p) = fill(f(p), 3) - - # Scalar integrand + # Integrand & Solution + integrand(p) = 1.0u"A" sol = Meshes.measure(triangle) * u"A" - @test integral(f, triangle, GaussLegendre(100)) ≈ sol - @test integral(f, triangle, GaussKronrod()) ≈ sol - @test integral(f, triangle, HAdaptiveCubature()) ≈ sol - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, triangle, GaussLegendre(100)) ≈ vsol - @test integral(fv, triangle, GaussKronrod()) ≈ vsol - @test integral(fv, triangle, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, triangle) - @test surfaceintegral(f, triangle) ≈ sol - @test_throws "not supported" volumeintegral(f, triangle) + # Package and run tests + testable = TestableGeometry(integrand, triangle, solution) + runtests(testable, SupportStatus(:surface)) end From d197f8378b9e1b75e0479017a93fd34bcaf40c28 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 11:51:27 -0500 Subject: [PATCH 17/47] Implement for Torus --- test/combinations.jl | 27 ++++++++------------------- 1 file changed, 8 insertions(+), 19 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index bd2d5f6b..4452f3ce 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -807,30 +807,19 @@ end @test volumeintegral(f, tetrahedron, HAdaptiveCubature()) ≈ sol end -@testitem "Meshes.Torus" setup=[Setup] begin +@testitem "Meshes.Torus" setup=[Combinations] begin + # Geometry origin = Point(0, 0, 0) ẑ = Vec(0, 0, 1) torus = Torus(origin, ẑ, 3.5, 1.25) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = Meshes.measure(torus) - @test integral(f, torus, GaussLegendre(100)) ≈ sol - @test integral(f, torus, GaussKronrod()) ≈ sol - @test integral(f, torus, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, torus, GaussLegendre(100)) ≈ vsol - @test integral(fv, torus, GaussKronrod()) ≈ vsol - @test integral(fv, torus, HAdaptiveCubature()) ≈ vsol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(torus) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, torus) - @test surfaceintegral(f, torus) ≈ sol - @test_throws "not supported" volumeintegral(f, torus) + # Package and run tests + testable = TestableGeometry(integrand, torus, solution) + runtests(testable, SupportStatus(:surface)) end @testitem "Meshes.Triangle" setup=[Combinations] begin From c58cc68d67c9cb93128facf0987f2c4f4cd9a67d Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 12:16:40 -0500 Subject: [PATCH 18/47] Implement Spheres --- test/combinations.jl | 85 ++++++++++++++------------------------------ 1 file changed, 26 insertions(+), 59 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 4452f3ce..8a1dfb95 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -716,95 +716,62 @@ end @test_throws "not supported" volumeintegral(f, segment) end -@testitem "Meshes.Sphere 2D" setup=[Setup] begin +@testitem "Meshes.Sphere 2D" setup=[Combinations] begin + # Geometry origin = Point(0, 0) radius = 4.4 sphere = Sphere(origin, radius) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(-r^2) + exp(-r^2) * u" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = 2π * radius * exp(-radius^2) * u"m" - @test integral(f, sphere, GaussLegendre(100)) ≈ sol - @test integral(f, sphere, GaussKronrod()) ≈ sol - @test integral(f, sphere, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, sphere, GaussLegendre(100)) ≈ vsol - @test integral(fv, sphere, GaussKronrod()) ≈ vsol - @test integral(fv, sphere, HAdaptiveCubature()) ≈ vsol + solution = 2π * radius * exp(-radius^2) * u"m" - # Integral aliases - @test lineintegral(f, sphere) ≈ sol - @test_throws "not supported" surfaceintegral(f, sphere) - @test_throws "not supported" volumeintegral(f, sphere) + # Package and run tests + testable = TestableGeometry(integrand, sphere, solution) + runtests(testable, SupportStatus(:line)) end -@testitem "Meshes.Sphere 3D" setup=[Setup] begin +@testitem "Meshes.Sphere 3D" setup=[Combinations] begin using CoordRefSystems: Cartesian, Spherical + # Geometry center = Point(1, 2, 3) radius = 4.4u"m" sphere = Sphere(center, radius) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} rθφ = convert(Spherical, Cartesian((p - center)...)) r = ustrip(rθφ.r) θ = ustrip(rθφ.θ) φ = ustrip(rθφ.ϕ) - sin(φ)^2 + cos(θ)^2 + (sin(φ)^2 + cos(θ)^2) * u"A" end - fv(p) = fill(f(p), 3) + solution = ((2π * radius^2) + (4π / 3 * radius^2)) * u"A" - # Scalar integrand - sol = (2π * radius^2) + (4π / 3 * radius^2) - @test integral(f, sphere, GaussLegendre(100)) ≈ sol - @test integral(f, sphere, GaussKronrod()) ≈ sol - @test integral(f, sphere, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, sphere, GaussLegendre(100)) ≈ vsol - @test integral(fv, sphere, GaussKronrod()) ≈ vsol - @test integral(fv, sphere, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, sphere) - @test surfaceintegral(f, sphere) ≈ sol - @test_throws "not supported" volumeintegral(f, sphere) + # Package and run tests + testable = TestableGeometry(integrand, sphere, solution) + runtests(testable, SupportStatus(:surface)) end -@testitem "Meshes.Tetrahedron" setup=[Setup] begin +@testitem "Meshes.Tetrahedron" setup=[Combinations] begin + # Geometry pt_n = Point(0, 3, 0) pt_w = Point(-7, 0, 0) pt_e = Point(8, 0, 0) ẑ = Vec(0, 0, 1) tetrahedron = Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ) - f(p) = 1.0u"A" - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = Meshes.measure(tetrahedron) * u"A" - @test integral(f, tetrahedron, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, tetrahedron, GaussKronrod()) - @test integral(f, tetrahedron, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, tetrahedron, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, tetrahedron, GaussKronrod()) - @test integral(fv, tetrahedron, HAdaptiveCubature()) ≈ vsol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(tetrahedron) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, tetrahedron) - @test_throws "not supported" surfaceintegral(f, tetrahedron) - @test volumeintegral(f, tetrahedron, HAdaptiveCubature()) ≈ sol + # Package and run tests + testable = TestableGeometry(integrand, tetrahedron, solution) + runtests(testable, SupportStatus(:volume)) end @testitem "Meshes.Torus" setup=[Combinations] begin From 936d38b93abe3bce8e37f7f4406b3408510ac842 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 12:54:36 -0500 Subject: [PATCH 19/47] Fix typos --- test/combinations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 8a1dfb95..d7df57f1 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -725,7 +725,7 @@ end # Integrand & Solution function integrand(p::P) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(-r^2) * u" + exp(-r^2) * u"A" end solution = 2π * radius * exp(-radius^2) * u"m" @@ -798,7 +798,7 @@ end # Integrand & Solution integrand(p) = 1.0u"A" - sol = Meshes.measure(triangle) * u"A" + solution = Meshes.measure(triangle) * u"A" # Package and run tests testable = TestableGeometry(integrand, triangle, solution) From e06af4edc076c5220b4aa22760145c0ed8095abe Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:04:07 -0500 Subject: [PATCH 20/47] Implement for Rope and Segment --- test/combinations.jl | 72 +++++++++++++++----------------------------- 1 file changed, 24 insertions(+), 48 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index d7df57f1..54543e8d 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -652,68 +652,44 @@ end @test_throws "not supported" volumeintegral(f, rope) end -@testitem "Meshes.Rope" setup=[Setup] begin - pt_a = Point(0.0u"m", 0.0u"m", 0.0u"m") - pt_b = Point(1.0u"m", 0.0u"m", 0.0u"m") - pt_c = Point(1.0u"m", 1.0u"m", 0.0u"m") - pt_d = Point(1.0u"m", 1.0u"m", 1.0u"m") +@testitem "Meshes.Rope" setup=[Combinations] begin + # Geometry + pt_a = Point(0, 0, 0) + pt_b = Point(1, 0, 0) + pt_c = Point(1, 1, 0) + pt_d = Point(1, 1, 1) rope = Rope(pt_a, pt_b, pt_c, pt_d) - function f(p::P) where {P <: Meshes.Point} - x, y, z = (p.coords.x, p.coords.y, p.coords.z) - (x + 2y + 3z) * u"A/m^2" + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} + x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) + (x + 2y + 3z) * u"A" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = 7.0u"A" - @test integral(f, rope, GaussLegendre(100)) ≈ sol - @test integral(f, rope, GaussKronrod()) ≈ sol - @test integral(f, rope, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, rope, GaussLegendre(100)) ≈ vsol - @test integral(fv, rope, GaussKronrod()) ≈ vsol - @test integral(fv, rope, HAdaptiveCubature()) ≈ vsol + solution = 7.0u"A*m" - # Integral aliases - @test lineintegral(f, rope) ≈ sol - @test_throws "not supported" surfaceintegral(f, rope) - @test_throws "not supported" volumeintegral(f, rope) + # Package and run tests + testable = TestableGeometry(integrand, rope, solution) + runtests(testable, SupportStatus(:line)) end -@testitem "Meshes.Segment" setup=[Setup] begin +@testitem "Meshes.Segment" setup=[Combinations] begin # Connect a line segment from the origin to an arbitrary point on the unit sphere φ, θ = (7pi / 6, pi / 3) # Arbitrary spherical angles - pt_a = Point(0.0u"m", 0.0u"m", 0.0u"m") - pt_b = Point(sin(θ) * cos(φ) * u"m", sin(θ) * sin(φ) * u"m", cos(θ) * u"m") + pt_a = Point(0, 0, 0) + pt_b = Point(sin(θ) * cos(φ), sin(θ) * sin(φ), cos(θ)) segment = Segment(pt_a, pt_b) + # Integrand & Solution a, b = (7.1, 4.6) # arbitrary constants > 0 - - function f(p::P) where {P <: Meshes.Point} + function integrand(p::P; a = a, b = b) where {P <: Meshes.Point} r = ustrip(u"m", norm(to(p))) - exp(r * log(a) + (1 - r) * log(b)) + exp(r * log(a) + (1 - r) * log(b)) * u"A" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = (a - b) / (log(a) - log(b)) * u"m" - @test integral(f, segment, GaussLegendre(100)) ≈ sol - @test integral(f, segment, GaussKronrod()) ≈ sol - @test integral(f, segment, HAdaptiveCubature()) ≈ sol - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, segment, GaussLegendre(100)) ≈ vsol - @test integral(fv, segment, GaussKronrod()) ≈ vsol - @test integral(fv, segment, HAdaptiveCubature()) ≈ vsol + solution = ((a - b) / (log(a) - log(b))) * u"A*m" - # Integral aliases - @test lineintegral(f, segment) ≈ sol - @test_throws "not supported" surfaceintegral(f, segment) - @test_throws "not supported" volumeintegral(f, segment) + # Package and run tests + testable = TestableGeometry(integrand, segment, solution) + runtests(testable, SupportStatus(:line)) end @testitem "Meshes.Sphere 2D" setup=[Combinations] begin From 3bd5dd85f40c7985856d4bccf9068c18d2ae3510 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:05:25 -0500 Subject: [PATCH 21/47] Fix units --- test/combinations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index 54543e8d..b8ee1f75 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -703,7 +703,7 @@ end r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end - solution = 2π * radius * exp(-radius^2) * u"m" + solution = 2π * radius * exp(-radius^2) * u"A*m" # Package and run tests testable = TestableGeometry(integrand, sphere, solution) From b1043b79ab0e3c895b09192a1d5974045bbe1d02 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:12:47 -0500 Subject: [PATCH 22/47] Implement for Ring --- test/combinations.jl | 51 +++++++++++++++++--------------------------- 1 file changed, 20 insertions(+), 31 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index b8ee1f75..fe45f4f9 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -621,44 +621,33 @@ end runtests(testable, SupportStatus(:line)) end -@testitem "Meshes.Ring" setup=[Setup] begin - pt_a = Point(0.0u"m", 0.0u"m", 0.0u"m") - pt_b = Point(1.0u"m", 0.0u"m", 0.0u"m") - pt_c = Point(1.0u"m", 1.0u"m", 0.0u"m") - pt_d = Point(1.0u"m", 1.0u"m", 1.0u"m") - rope = Ring(pt_a, pt_b, pt_c, pt_d, pt_c, pt_b) +@testitem "Meshes.Ring" setup=[Combinations] begin + # Geometry + a = Point(0, 0, 0) + b = Point(1, 0, 0) + c = Point(1, 1, 0) + d = Point(1, 1, 1) + ring = Ring(a, b, c, d, c, b) - function f(p::P) where {P <: Meshes.Point} - x, y, z = (p.coords.x, p.coords.y, p.coords.z) - (x + 2y + 3z) * u"A/m^2" + # Integrand & Solution + function integrand(p::P) where {P <: Meshes.Point} + x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) + (x + 2y + 3z) * u"A" end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = 14.0u"A" - @test integral(f, rope, GaussLegendre(100)) ≈ sol - @test integral(f, rope, GaussKronrod()) ≈ sol - @test integral(f, rope, HAdaptiveCubature()) ≈ sol + solution = 14.0u"A*m" - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, rope, GaussLegendre(100)) ≈ vsol - @test integral(fv, rope, GaussKronrod()) ≈ vsol - @test integral(fv, rope, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test lineintegral(f, rope) ≈ sol - @test_throws "not supported" surfaceintegral(f, rope) - @test_throws "not supported" volumeintegral(f, rope) + # Package and run tests + testable = TestableGeometry(integrand, ring, solution) + runtests(testable, SupportStatus(:line)) end @testitem "Meshes.Rope" setup=[Combinations] begin # Geometry - pt_a = Point(0, 0, 0) - pt_b = Point(1, 0, 0) - pt_c = Point(1, 1, 0) - pt_d = Point(1, 1, 1) - rope = Rope(pt_a, pt_b, pt_c, pt_d) + a = Point(0, 0, 0) + b = Point(1, 0, 0) + c = Point(1, 1, 0) + d = Point(1, 1, 1) + rope = Rope(a, b, c, d) # Integrand & Solution function integrand(p::P) where {P <: Meshes.Point} From 78a65d03fefd06634ba31f44b23e8a56abf0c3fc Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:14:33 -0500 Subject: [PATCH 23/47] Absorb integrand parameters --- test/combinations.jl | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index fe45f4f9..5428fab9 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -103,7 +103,7 @@ end #testsnippet ball = Ball(origin, radius) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) exp(-r^2) end @@ -123,7 +123,7 @@ end ball = Ball(center, radius) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) offset = p - center ur = norm(offset) r = ustrip(u"m", ur) @@ -143,7 +143,7 @@ end curve = BezierCurve([Point(t, sin(t), 0) for t in range(-π, π, length = 361)]) # Integrand - function f(p::P) where {P <: Meshes.Point} + function f(p::Meshes.Point) ux = ustrip(p.coords.x) (1 / sqrt(1 + cos(ux)^2)) * u"Ω/m" end @@ -205,7 +205,7 @@ end a = π box = Box(Point(0, 0), Point(a, a)) - function f(p::P) where {P <: Meshes.Point} + function f(p::Meshes.Point) x, y = ustrip.((p.coords.x, p.coords.y)) (sqrt(a^2 - x^2) + sqrt(a^2 - y^2)) * u"Ω/m^2" end @@ -236,7 +236,7 @@ end a = π box = Box(Point(0, 0, 0), Point(a, a, a)) - function f(p::P) where {P <: Meshes.Point} + function f(p::Meshes.Point) x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3" end @@ -264,7 +264,7 @@ end a = π box = Box(Point(0, 0, 0, 0), Point(a, a, a, a)) - function f(p::P) where {P <: Meshes.Point} + function f(p::Meshes.Point) x1, x2, x3, x4 = ustrip.(to(p).coords) σ(x) = sqrt(a^2 - x^2) (σ(x1) + σ(x2) + σ(x3) + σ(x4)) * u"Ω/m^4" @@ -298,7 +298,7 @@ end circle = Circle(plane, radius) # Integrand - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) offset = p - center r = ustrip(u"m", norm(offset)) exp(-r^2) @@ -404,7 +404,7 @@ end disk = Disk(plane, radius) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) offset = p - center r = ustrip(u"m", norm(offset)) exp(-r^2) * u"A" @@ -514,7 +514,7 @@ end line = Line(a, b) # Integrand & solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end @@ -552,7 +552,7 @@ end curve_polar = ParametrizedCurve(t -> Point(Polar(radius, t)), (0.0, 2π)) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) ur = norm(to(p)) r = ustrip(u"m", ur) exp(-r^2) * u"A" @@ -574,7 +574,7 @@ end plane = Plane(p, v) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end @@ -592,7 +592,7 @@ end quadrangle = Quadrangle((-1.0, 0.0), (-1.0, 1.0), (1.0, 1.0), (1.0, 0.0)) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end @@ -610,7 +610,7 @@ end ray = Ray(a, v) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end @@ -630,7 +630,7 @@ end ring = Ring(a, b, c, d, c, b) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) (x + 2y + 3z) * u"A" end @@ -650,7 +650,7 @@ end rope = Rope(a, b, c, d) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) (x + 2y + 3z) * u"A" end @@ -688,7 +688,7 @@ end sphere = Sphere(origin, radius) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) exp(-r^2) * u"A" end @@ -708,7 +708,7 @@ end sphere = Sphere(center, radius) # Integrand & Solution - function integrand(p::P) where {P <: Meshes.Point} + function integrand(p::Meshes.Point) rθφ = convert(Spherical, Cartesian((p - center)...)) r = ustrip(rθφ.r) θ = ustrip(rθφ.θ) From 0f1612173cc7b9d318051bc42f14a8850c61a288 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:21:00 -0500 Subject: [PATCH 24/47] Code re-org and remove Box-4D --- test/combinations.jl | 78 ++++++++++++++------------------------------ 1 file changed, 24 insertions(+), 54 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 5428fab9..ec880186 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -120,6 +120,7 @@ end # Geometry center = Point(1, 2, 3) radius = 2.8u"m" + r = ustrip(u"m", radius) ball = Ball(center, radius) # Integrand & Solution @@ -129,9 +130,7 @@ end r = ustrip(u"m", ur) exp(-r^2) end - solution = let r = ustrip(u"m", radius) - (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" - end + solution = (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" # Package and run tests testable = TestableGeometry(integrand, ball, solution) @@ -145,12 +144,11 @@ end # Integrand function f(p::Meshes.Point) ux = ustrip(p.coords.x) - (1 / sqrt(1 + cos(ux)^2)) * u"Ω/m" + (1 / sqrt(1 + cos(ux)^2)) * u"Ω" end - fv(p) = fill(f(p), 3) + solution = 2π * u"Ω*m" # Scalar integrand - sol = 2π * u"Ω" @test integral(f, curve, GaussLegendre(100))≈sol rtol=0.5e-2 @test integral(f, curve, GaussKronrod())≈sol rtol=0.5e-2 @test integral(f, curve, HAdaptiveCubature())≈sol rtol=0.5e-2 @@ -171,34 +169,33 @@ end end @testitem "Meshes.Box 1D" setup=[Setup] begin + # Geometry a = π box = Box(Point(0), Point(a)) - struct Integrand end - fc = Integrand() - - function (::Integrand)(p::Meshes.Point) + # Integrand & Solution + function f(p::Meshes.Point) t = ustrip(p.coords.x) - sqrt(a^2 - t^2) * u"Ω/m" + sqrt(a^2 - t^2) * u"Ω" end - fcv(p) = fill(fc(p), 3) + fv(p) = fill(f(p), 3) + solution = π * a^2 / 4 * u"Ω*m" + vsol = fill(solution, 3) # Scalar integrand - sol = π * a^2 / 4 * u"Ω" - @test integral(fc, box, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(fc, box, GaussKronrod()) ≈ sol - @test integral(fc, box, HAdaptiveCubature()) ≈ sol + @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 + @test integral(f, box, GaussKronrod()) ≈ sol + @test integral(f, box, HAdaptiveCubature()) ≈ sol # Vector integrand - vsol = fill(sol, 3) - @test integral(fcv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fcv, box, GaussKronrod()) ≈ vsol - @test integral(fcv, box, HAdaptiveCubature()) ≈ vsol + @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 + @test integral(fv, box, GaussKronrod()) ≈ vsol + @test integral(fv, box, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(fc, box) ≈ sol - @test_throws "not supported" surfaceintegral(fc, box) - @test_throws "not supported" volumeintegral(fc, box) + @test lineintegral(f, box) ≈ sol + @test_throws "not supported" surfaceintegral(f, box) + @test_throws "not supported" volumeintegral(f, box) end @testitem "Meshes.Box 2D" setup=[Setup] begin @@ -233,23 +230,25 @@ end end @testitem "Meshes.Box 3D" setup=[Setup] begin + # Geometry a = π box = Box(Point(0, 0, 0), Point(a, a, a)) + # Integrand & Solution function f(p::Meshes.Point) x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3" end fv(p) = fill(f(p), 3) + sol = 3a^2 * (π * a^2 / 4) * u"Ω" + vsol = fill(sol, 3) # Scalar integrand - sol = 3a^2 * (π * a^2 / 4) * u"Ω" @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 @test_throws "not supported" integral(f, box, GaussKronrod()) @test integral(f, box, HAdaptiveCubature()) ≈ sol # Vector integrand - vsol = fill(sol, 3) @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 @test_throws "not supported" integral(fv, box, GaussKronrod()) @test integral(fv, box, HAdaptiveCubature()) ≈ vsol @@ -260,35 +259,6 @@ end @test volumeintegral(f, box) ≈ sol end -@testitem "Meshes.Box 4D" tags=[:extended] setup=[Setup] begin - a = π - box = Box(Point(0, 0, 0, 0), Point(a, a, a, a)) - - function f(p::Meshes.Point) - x1, x2, x3, x4 = ustrip.(to(p).coords) - σ(x) = sqrt(a^2 - x^2) - (σ(x1) + σ(x2) + σ(x3) + σ(x4)) * u"Ω/m^4" - end - fv(p) = fill(f(p), 3) - - # Scalar integrand - sol = 4a^3 * (π * a^2 / 4) * u"Ω" - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test_throws "not supported" integral(f, box, GaussKronrod()) - @test integral(f, box, HAdaptiveCubature(rtol = 1e-6))≈sol rtol=1e-6 - - # Vector integrand - vsol = fill(sol, 3) - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test_throws "not supported" integral(fv, box, GaussKronrod()) - @test integral(fv, box, HAdaptiveCubature(rtol = 1e-6))≈vsol rtol=1e-6 - - # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test_throws "not supported" surfaceintegral(f, box) - @test_throws "not supported" volumeintegral(f, box) -end - @testitem "Meshes.Circle" setup=[Combinations] begin # Geometry center = Point(1, 2, 3) From d65730dbd8a05a300e27bea653011876c17defa4 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:27:07 -0500 Subject: [PATCH 25/47] Code cleanup and add Unitful integrands --- test/combinations.jl | 51 ++++++++++++++++++++++---------------------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index ec880186..c83a9340 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -105,9 +105,9 @@ end #testsnippet # Integrand & Solution function integrand(p::Meshes.Point) r = ustrip(u"m", norm(to(p))) - exp(-r^2) + exp(-r^2) * u"A" end - solution = (π - π * exp(-radius^2)) * u"m^2" + solution = (π - π * exp(-radius^2)) * u"A*m^2" # Package and run tests testable = TestableGeometry(integrand, ball, solution) @@ -128,9 +128,9 @@ end offset = p - center ur = norm(offset) r = ustrip(u"m", ur) - exp(-r^2) + exp(-r^2) * u"A" end - solution = (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" + solution = (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"A*m^3" # Package and run tests testable = TestableGeometry(integrand, ball, solution) @@ -267,15 +267,13 @@ end radius = 4.4 circle = Circle(plane, radius) - # Integrand + # Integrand & Solution function integrand(p::Meshes.Point) offset = p - center r = ustrip(u"m", norm(offset)) - exp(-r^2) + exp(-r^2) * u"A" end - - # Scalar integrand - solution = 2π * radius * exp(-radius^2) * u"m" + solution = 2π * radius * exp(-radius^2) * u"A*m" # Package and run tests testable = TestableGeometry(integrand, circle, solution) @@ -292,10 +290,8 @@ end apex = Point(0.0u"m", 0.0u"m", h) cone = Cone(base, apex) - # Integrand + # Integrand & Solution integrand(p) = 1.0u"A" - - # Solution solution = (π * r^2 * h / 3) * u"A" # Package and run tests @@ -304,6 +300,7 @@ end end @testitem "Meshes.ConeSurface" setup=[Setup] begin + # Geometry r = 2.5u"m" h = 2.5u"m" origin = Point(0, 0, 0) @@ -312,19 +309,18 @@ end apex = Point(0.0u"m", 0.0u"m", h) cone = ConeSurface(base, apex) - f(p) = 1.0 + # Integrand & Solution + f(p) = 1.0u"A" fv(p) = fill(f(p), 3) - - _area_cone_rightcircular(h, r) = (π * r^2) + (π * r * hypot(h, r)) + sol = ((π * r^2) + (π * r * hypot(h, r))) * u"A" + vsol = fill(sol, 3) # Scalar integrand - sol = _area_cone_rightcircular(h, r) @test integral(f, cone, GaussLegendre(100))≈sol rtol=1e-6 @test integral(f, cone, GaussKronrod())≈sol rtol=1e-6 @test integral(f, cone, HAdaptiveCubature()) ≈ sol # Vector integrand - vsol = fill(sol, 3) @test integral(fv, cone, GaussLegendre(100))≈vsol rtol=1e-6 @test integral(fv, cone, GaussKronrod())≈vsol rtol=1e-6 @test integral(fv, cone, HAdaptiveCubature()) ≈ vsol @@ -387,22 +383,24 @@ end end @testitem "Meshes.Ellipsoid" setup=[Setup] begin + # Geometry origin = Point(0, 0, 0) radii = (1.0, 2.0, 0.5) ellipsoid = Ellipsoid(radii, origin) - f(p) = 1.0 + # Integrand & Solution + f(p) = 1.0u"A" fv(p) = fill(f(p), 3) + sol = Meshes.measure(ellipsoid) * u"A" + vsol = fill(sol, 3) # Tolerances are higher due to `measure` being only an approximation # Scalar integrand - sol = Meshes.measure(ellipsoid) @test integral(f, ellipsoid, GaussLegendre(100))≈sol rtol=1e-2 @test integral(f, ellipsoid, GaussKronrod())≈sol rtol=1e-2 @test integral(f, ellipsoid, HAdaptiveCubature())≈sol rtol=1e-2 # Vector integrand - vsol = fill(sol, 3) @test integral(fv, ellipsoid, GaussLegendre(100))≈vsol rtol=1e-2 @test integral(fv, ellipsoid, GaussKronrod())≈vsol rtol=1e-2 @test integral(fv, ellipsoid, HAdaptiveCubature())≈vsol rtol=1e-2 @@ -428,25 +426,26 @@ end disk_top = Disk(plane_top, r_top) frustum = FrustumSurface(disk_bot, disk_top) - f(p) = 1.0 + # Integrand & Solution + f(p) = 1.0u"A" fv(p) = fill(f(p), 3) - _area_base(r) = π * r^2 _area_cone_walls(h, r) = π * r * hypot(h, r) - - # Scalar integrand sol = let area_walls_projected = _area_cone_walls(cone_h, r_bot) area_walls_missing = _area_cone_walls(0.5cone_h, r_top) area_walls = area_walls_projected - area_walls_missing - area_walls + _area_base(r_top) + _area_base(r_bot) + area_total = area_walls + _area_base(r_top) + _area_base(r_bot) + area_total * u"A" end + vsol = fill(sol, 3) + + # Scalar integrand @test integral(f, frustum, GaussLegendre(100))≈sol rtol=1e-6 @test integral(f, frustum, GaussKronrod())≈sol rtol=1e-6 @test integral(f, frustum, HAdaptiveCubature()) ≈ sol # Vector integrand - vsol = fill(sol, 3) @test integral(fv, frustum, GaussLegendre(100))≈vsol rtol=1e-6 @test integral(fv, frustum, GaussKronrod())≈vsol rtol=1e-6 @test integral(fv, frustum, HAdaptiveCubature()) ≈ vsol From 43edac000baf06bf6ecb100c4c7aa51d09131f9b Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 13:29:11 -0500 Subject: [PATCH 26/47] Code re-org in Hexahedron --- test/combinations.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index c83a9340..6b44d109 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -452,20 +452,22 @@ end end @testitem "Meshes.Hexahedron" setup=[Setup] begin + # Geometry hexahedron = Hexahedron(Point(0, 0, 0), Point(2, 0, 0), Point(2, 2, 0), Point(0, 2, 0), Point(0, 0, 2), Point(1, 0, 2), Point(1, 1, 2), Point(0, 1, 2)) - f(p) = 1.0 + # Integrand & Solution + f(p) = 1.0u"A" fv(p) = fill(f(p), 3) + sol = Meshes.measure(hexahedron) * u"A" + vsol = fill(sol, 3) # Scalar integrand - sol = Meshes.measure(hexahedron) @test integral(f, hexahedron, GaussLegendre(100)) ≈ sol @test_throws "not supported" integral(f, hexahedron, GaussKronrod())≈sol @test integral(f, hexahedron, HAdaptiveCubature()) ≈ sol # Vector integrand - vsol = fill(sol, 3) @test integral(fv, hexahedron, GaussLegendre(100)) ≈ vsol @test_throws "not supported" integral(fv, hexahedron, GaussKronrod())≈vsol @test integral(fv, hexahedron, HAdaptiveCubature()) ≈ vsol From b8decd59985409dd15673672b7684f0d28b67c88 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 13:40:40 -0500 Subject: [PATCH 27/47] Typo fix --- test/combinations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index c83a9340..26f4fdcf 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -146,7 +146,8 @@ end ux = ustrip(p.coords.x) (1 / sqrt(1 + cos(ux)^2)) * u"Ω" end - solution = 2π * u"Ω*m" + sol = 2π * u"Ω*m" + vsol = fill(sol, 3) # Scalar integrand @test integral(f, curve, GaussLegendre(100))≈sol rtol=0.5e-2 @@ -154,7 +155,6 @@ end @test integral(f, curve, HAdaptiveCubature())≈sol rtol=0.5e-2 # Vector integrand - vsol = fill(sol, 3) @test integral(fv, curve, GaussLegendre(100))≈vsol rtol=0.5e-2 @test integral(fv, curve, GaussKronrod())≈vsol rtol=0.5e-2 @test integral(fv, curve, HAdaptiveCubature())≈vsol rtol=0.5e-2 From 452278b4afb71453a00502d4e8be21decd30c2d1 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 14:04:54 -0500 Subject: [PATCH 28/47] Bugfixes --- test/combinations.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 26f4fdcf..4649e437 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -146,6 +146,7 @@ end ux = ustrip(p.coords.x) (1 / sqrt(1 + cos(ux)^2)) * u"Ω" end + fv(p) = fill(f(p), 3) sol = 2π * u"Ω*m" vsol = fill(sol, 3) @@ -179,8 +180,8 @@ end sqrt(a^2 - t^2) * u"Ω" end fv(p) = fill(f(p), 3) - solution = π * a^2 / 4 * u"Ω*m" - vsol = fill(solution, 3) + sol = π * a^2 / 4 * u"Ω*m" + vsol = fill(sol, 3) # Scalar integrand @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 From 2d6a6ea722e0666a7c29bc1f0ee8836679d78cb4 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 14:18:24 -0500 Subject: [PATCH 29/47] Add rtol kwarg to runtests --- test/combinations.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 98840763..143a89c0 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -53,12 +53,12 @@ end end - function runtests(testable::TestableGeometry, supports::SupportStatus) + function runtests(testable::TestableGeometry, supports::SupportStatus; rtol=sqrt(eps())) # Test alias functions for alias in (lineintegral, surfaceintegral, volumeintegral) alias_symbol = first(methods(alias)).name if getfield(supports, alias_symbol) - @test alias(testable.integrand, testable.geometry) ≈ testable.solution + @test alias(testable.integrand, testable.geometry) ≈ testable.solution rtol=rtol else @test_throws "not supported" alias(testable.integrand, testable.geometry) end @@ -75,16 +75,16 @@ if supported # Scalar integrand sol = testable.solution - @test integral(testable.integrand, testable.geometry, rule) ≈ sol + @test integral(testable.integrand, testable.geometry, rule) ≈ sol rtol=rtol # Callable integrand f = Callable(testable.integrand) - @test integral(f, testable.geometry, rule) ≈ sol + @test integral(f, testable.geometry, rule) ≈ sol rtol=rtol # Vector integrand fv(p) = fill(testable.integrand(p), 3) sol_v = fill(testable.solution, 3) - @test integral(fv, testable.geometry, rule) ≈ sol_v + @test integral(fv, testable.geometry, rule) ≈ sol_v rtol=rtol else @test_throws "not supported" integral(testable.integrand, testable.geometry, rule) end # if From b0529af0358d848c3ecaec99b3ead1a49e28b9fe Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 14:28:44 -0500 Subject: [PATCH 30/47] Implement for BezierCurve and Box-1D --- test/combinations.jl | 51 +++++++++++--------------------------------- 1 file changed, 13 insertions(+), 38 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 143a89c0..33fba8bb 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -137,66 +137,41 @@ end runtests(testable, SupportStatus(:volume)) end -@testitem "Meshes.BezierCurve" setup=[Setup] begin +@testitem "Meshes.BezierCurve" setup=[Combinations] begin # Geometry curve = BezierCurve([Point(t, sin(t), 0) for t in range(-π, π, length = 361)]) # Integrand - function f(p::Meshes.Point) + function integrand(p::Meshes.Point) ux = ustrip(p.coords.x) (1 / sqrt(1 + cos(ux)^2)) * u"Ω" end - fv(p) = fill(f(p), 3) - sol = 2π * u"Ω*m" - vsol = fill(sol, 3) - - # Scalar integrand - @test integral(f, curve, GaussLegendre(100))≈sol rtol=0.5e-2 - @test integral(f, curve, GaussKronrod())≈sol rtol=0.5e-2 - @test integral(f, curve, HAdaptiveCubature())≈sol rtol=0.5e-2 - - # Vector integrand - @test integral(fv, curve, GaussLegendre(100))≈vsol rtol=0.5e-2 - @test integral(fv, curve, GaussKronrod())≈vsol rtol=0.5e-2 - @test integral(fv, curve, HAdaptiveCubature())≈vsol rtol=0.5e-2 + solution = 2π * u"Ω*m" - # Integral aliases - @test lineintegral(f, curve)≈sol rtol=0.5e-2 - @test_throws "not supported" surfaceintegral(f, curve) - @test_throws "not supported" volumeintegral(f, curve) + # Package and run tests + testable = TestableGeometry(integrand, curve, solution) + runtests(testable, SupportStatus(:line); rtol=0.5e-2) + # TODO move this to Differentiation testitem # Check Bezier-specific jacobian bounds @test_throws DomainError jacobian(curve, [1.1]) end -@testitem "Meshes.Box 1D" setup=[Setup] begin +@testitem "Meshes.Box 1D" setup=[Combinations] begin # Geometry a = π box = Box(Point(0), Point(a)) # Integrand & Solution - function f(p::Meshes.Point) + function integrand(p::Meshes.Point) t = ustrip(p.coords.x) sqrt(a^2 - t^2) * u"Ω" end - fv(p) = fill(f(p), 3) - sol = π * a^2 / 4 * u"Ω*m" - vsol = fill(sol, 3) + solution = π * a^2 / 4 * u"Ω*m" - # Scalar integrand - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(f, box, GaussKronrod()) ≈ sol - @test integral(f, box, HAdaptiveCubature()) ≈ sol - - # Vector integrand - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, box, GaussKronrod()) ≈ vsol - @test integral(fv, box, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test lineintegral(f, box) ≈ sol - @test_throws "not supported" surfaceintegral(f, box) - @test_throws "not supported" volumeintegral(f, box) + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable, SupportStatus(:line); rtol=1e-6) end @testitem "Meshes.Box 2D" setup=[Setup] begin From 40d67dbb83dcdb9e569d5906ad2e765bc07d6f5f Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 14:39:05 -0500 Subject: [PATCH 31/47] Implement for Box 2D and 3D --- test/combinations.jl | 52 ++++++++++++-------------------------------- 1 file changed, 14 insertions(+), 38 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 33fba8bb..ae82cc37 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -174,65 +174,41 @@ end runtests(testable, SupportStatus(:line); rtol=1e-6) end -@testitem "Meshes.Box 2D" setup=[Setup] begin +@testitem "Meshes.Box 2D" setup=[Combinations] begin a = π box = Box(Point(0, 0), Point(a, a)) - function f(p::Meshes.Point) + # Integrand & Solution + function integrand(p::Meshes.Point) x, y = ustrip.((p.coords.x, p.coords.y)) (sqrt(a^2 - x^2) + sqrt(a^2 - y^2)) * u"Ω/m^2" end - fv(p) = fill(f(p), 3) + solution = 2a * (π * a^2 / 4) * u"Ω" - # Scalar integrand - sol = 2a * (π * 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 - - # 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test surfaceintegral(f, box) ≈ sol - @test_throws "not supported" volumeintegral(f, box) + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable, SupportStatus(:surface); rtol=1e-6) + # TODO move this to Differentiation testitem # Test jacobian with wrong number of parametric coordinates @test_throws ArgumentError jacobian(box, zeros(3)) end -@testitem "Meshes.Box 3D" setup=[Setup] begin +@testitem "Meshes.Box 3D" setup=[Combinations] begin # Geometry a = π box = Box(Point(0, 0, 0), Point(a, a, a)) # Integrand & Solution - function f(p::Meshes.Point) + function integrand(p::Meshes.Point) x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3" end - fv(p) = fill(f(p), 3) - sol = 3a^2 * (π * a^2 / 4) * u"Ω" - vsol = fill(sol, 3) - - # Scalar integrand - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test_throws "not supported" integral(f, box, GaussKronrod()) - @test integral(f, box, HAdaptiveCubature()) ≈ sol + solution = 3a^2 * (π * a^2 / 4) * u"Ω" - # Vector integrand - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test_throws "not supported" integral(fv, box, GaussKronrod()) - @test integral(fv, box, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test_throws "not supported" surfaceintegral(f, box) - @test volumeintegral(f, box) ≈ sol + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable, SupportStatus(:volume); rtol=1e-6) end @testitem "Meshes.Circle" setup=[Combinations] begin From 122542d47c95e1919acd1d4d3c6aff6060f954e5 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 14:44:42 -0500 Subject: [PATCH 32/47] Implement ConeSurface and Ellipsoid --- test/combinations.jl | 52 +++++++++++--------------------------------- 1 file changed, 13 insertions(+), 39 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index ae82cc37..c5620afc 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -251,10 +251,10 @@ end runtests(testable, SupportStatus(:volume)) end -@testitem "Meshes.ConeSurface" setup=[Setup] begin +@testitem "Meshes.ConeSurface" setup=[Combinations] begin # Geometry r = 2.5u"m" - h = 2.5u"m" + h = 3.5u"m" origin = Point(0, 0, 0) xy_plane = Plane(origin, Vec(0, 0, 1)) base = Disk(xy_plane, r) @@ -262,25 +262,12 @@ end cone = ConeSurface(base, apex) # Integrand & Solution - f(p) = 1.0u"A" - fv(p) = fill(f(p), 3) - sol = ((π * r^2) + (π * r * hypot(h, r))) * u"A" - vsol = fill(sol, 3) - - # Scalar integrand - @test integral(f, cone, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(f, cone, GaussKronrod())≈sol rtol=1e-6 - @test integral(f, cone, HAdaptiveCubature()) ≈ sol - - # Vector integrand - @test integral(fv, cone, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, cone, GaussKronrod())≈vsol rtol=1e-6 - @test integral(fv, cone, HAdaptiveCubature()) ≈ vsol + integrand(p) = 1.0u"A" + solution = ((π * r^2) + (π * r * hypot(h, r))) * u"A" - # Integral aliases - @test_throws "not supported" lineintegral(f, cone) - @test surfaceintegral(f, cone) ≈ sol - @test_throws "not supported" volumeintegral(f, cone) + # Package and run tests + testable = TestableGeometry(integrand, cone, solution) + runtests(testable, SupportStatus(:surface); rtol=1e-6) end @testitem "Meshes.Cylinder" setup=[Combinations] begin @@ -334,33 +321,20 @@ end runtests(testable, SupportStatus(:surface)) end -@testitem "Meshes.Ellipsoid" setup=[Setup] begin +@testitem "Meshes.Ellipsoid" setup=[Combinations] begin # Geometry origin = Point(0, 0, 0) radii = (1.0, 2.0, 0.5) ellipsoid = Ellipsoid(radii, origin) # Integrand & Solution - f(p) = 1.0u"A" - fv(p) = fill(f(p), 3) - sol = Meshes.measure(ellipsoid) * u"A" - vsol = fill(sol, 3) + integrand(p) = 1.0u"A" + solution = Meshes.measure(ellipsoid) * u"A" + # Package and run tests # Tolerances are higher due to `measure` being only an approximation - # Scalar integrand - @test integral(f, ellipsoid, GaussLegendre(100))≈sol rtol=1e-2 - @test integral(f, ellipsoid, GaussKronrod())≈sol rtol=1e-2 - @test integral(f, ellipsoid, HAdaptiveCubature())≈sol rtol=1e-2 - - # Vector integrand - @test integral(fv, ellipsoid, GaussLegendre(100))≈vsol rtol=1e-2 - @test integral(fv, ellipsoid, GaussKronrod())≈vsol rtol=1e-2 - @test integral(fv, ellipsoid, HAdaptiveCubature())≈vsol rtol=1e-2 - - # Integral aliases - @test_throws "not supported" lineintegral(f, ellipsoid) - @test surfaceintegral(f, ellipsoid)≈sol rtol=1e-2 - @test_throws "not supported" volumeintegral(f, ellipsoid) + testable = TestableGeometry(integrand, ellipsoid, solution) + runtests(testable, SupportStatus(:surface); rtol=1e-2) end @testitem "Meshes.FrustumSurface" setup=[Setup] begin From fefd49068561706aecc9d4927e554d34a60bc128 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 14:48:51 -0500 Subject: [PATCH 33/47] Implement for FrustumSurface and Hexahedron --- test/combinations.jl | 51 +++++++++++++------------------------------- 1 file changed, 15 insertions(+), 36 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index c5620afc..ae2fd299 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -337,9 +337,9 @@ end runtests(testable, SupportStatus(:surface); rtol=1e-2) end -@testitem "Meshes.FrustumSurface" setup=[Setup] begin - # Create a frustum whose radius halves at the top, - # i.e. the bottom half of a cone by height +@testitem "Meshes.FrustumSurface" setup=[Combinations] begin + # Geometry + # Create a frustum whose radius halves at the top, i.e. the bottom half of a cone r_bot = 2.5u"m" r_top = 1.25u"m" cone_h = 2π * u"m" @@ -353,55 +353,34 @@ end frustum = FrustumSurface(disk_bot, disk_top) # Integrand & Solution - f(p) = 1.0u"A" - fv(p) = fill(f(p), 3) + integrand(p) = 1.0u"A" _area_base(r) = π * r^2 _area_cone_walls(h, r) = π * r * hypot(h, r) - sol = let + solution = let area_walls_projected = _area_cone_walls(cone_h, r_bot) area_walls_missing = _area_cone_walls(0.5cone_h, r_top) area_walls = area_walls_projected - area_walls_missing area_total = area_walls + _area_base(r_top) + _area_base(r_bot) area_total * u"A" end - vsol = fill(sol, 3) - - # Scalar integrand - @test integral(f, frustum, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(f, frustum, GaussKronrod())≈sol rtol=1e-6 - @test integral(f, frustum, HAdaptiveCubature()) ≈ sol - # Vector integrand - @test integral(fv, frustum, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, frustum, GaussKronrod())≈vsol rtol=1e-6 - @test integral(fv, frustum, HAdaptiveCubature()) ≈ vsol + # Package and run tests + testable = TestableGeometry(integrand, frustum, solution) + runtests(testable, SupportStatus(:surface); rtol=1e-6) end -@testitem "Meshes.Hexahedron" setup=[Setup] begin +@testitem "Meshes.Hexahedron" setup=[Combinations] begin # Geometry hexahedron = Hexahedron(Point(0, 0, 0), Point(2, 0, 0), Point(2, 2, 0), Point(0, 2, 0), Point(0, 0, 2), Point(1, 0, 2), Point(1, 1, 2), Point(0, 1, 2)) # Integrand & Solution - f(p) = 1.0u"A" - fv(p) = fill(f(p), 3) - sol = Meshes.measure(hexahedron) * u"A" - vsol = fill(sol, 3) - - # Scalar integrand - @test integral(f, hexahedron, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, hexahedron, GaussKronrod())≈sol - @test integral(f, hexahedron, HAdaptiveCubature()) ≈ sol - - # Vector integrand - @test integral(fv, hexahedron, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, hexahedron, GaussKronrod())≈vsol - @test integral(fv, hexahedron, HAdaptiveCubature()) ≈ vsol - - # Integral aliases - @test_throws "not supported" lineintegral(f, hexahedron) - @test_throws "not supported" surfaceintegral(f, hexahedron) - @test volumeintegral(f, hexahedron) ≈ sol + integrand(p) = 1.0u"A" + solution = Meshes.measure(hexahedron) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, hexahedron, solution) + runtests(testable, SupportStatus(:volume)) end @testitem "Meshes.Line" setup=[Combinations] begin From 438193e1515e97bf3d7b222f81d3ed3239002d46 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 15:01:39 -0500 Subject: [PATCH 34/47] Move jacobian tests --- test/combinations.jl | 8 -------- test/utils.jl | 6 +++++- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index ae2fd299..4e015114 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -151,10 +151,6 @@ end # Package and run tests testable = TestableGeometry(integrand, curve, solution) runtests(testable, SupportStatus(:line); rtol=0.5e-2) - - # TODO move this to Differentiation testitem - # Check Bezier-specific jacobian bounds - @test_throws DomainError jacobian(curve, [1.1]) end @testitem "Meshes.Box 1D" setup=[Combinations] begin @@ -188,10 +184,6 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) runtests(testable, SupportStatus(:surface); rtol=1e-6) - - # TODO move this to Differentiation testitem - # Test jacobian with wrong number of parametric coordinates - @test_throws ArgumentError jacobian(box, zeros(3)) end @testitem "Meshes.Box 3D" setup=[Combinations] begin diff --git a/test/utils.jl b/test/utils.jl index 50cb4ea1..87ef390f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -19,7 +19,7 @@ @test _ones(Float32, 2) == (1.0f0, 1.0f0) end -@testitem "DifferentiationMethod" setup=[Setup] begin +@testitem "Differentiation" setup=[Setup] begin using MeshIntegrals: _default_diff_method # _default_diff_method @@ -29,6 +29,10 @@ end # FiniteDifference @test FiniteDifference().ε ≈ 1e-6 + + # Test jacobian with wrong number of parametric coordinates + box = Box(Point(0, 0), Point(a, a)) + @test_throws ArgumentError jacobian(box, zeros(3)) end @testitem "_ParametricGeometry" setup=[Setup] begin From d8d6a870f07e7558a65225d5aef6760139ef0183 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Thu, 5 Dec 2024 15:13:52 -0500 Subject: [PATCH 35/47] Bugfix --- test/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utils.jl b/test/utils.jl index 87ef390f..448d38de 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -31,7 +31,7 @@ end @test FiniteDifference().ε ≈ 1e-6 # Test jacobian with wrong number of parametric coordinates - box = Box(Point(0, 0), Point(a, a)) + box = Box(Point(0, 0), Point(1, 1)) @test_throws ArgumentError jacobian(box, zeros(3)) end From ebed2395c837455cbb71b7a9adc83041bee4f71a Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 5 Dec 2024 15:54:26 -0500 Subject: [PATCH 36/47] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/combinations.jl | 28 +++++++++++++++------------- test/utils.jl | 2 +- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 4e015114..5ed5bc95 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -53,12 +53,13 @@ end end - function runtests(testable::TestableGeometry, supports::SupportStatus; rtol=sqrt(eps())) + function runtests( + testable::TestableGeometry, supports::SupportStatus; rtol = sqrt(eps())) # Test alias functions for alias in (lineintegral, surfaceintegral, volumeintegral) alias_symbol = first(methods(alias)).name if getfield(supports, alias_symbol) - @test alias(testable.integrand, testable.geometry) ≈ testable.solution rtol=rtol + @test alias(testable.integrand, testable.geometry)≈testable.solution rtol=rtol else @test_throws "not supported" alias(testable.integrand, testable.geometry) end @@ -75,18 +76,19 @@ if supported # Scalar integrand sol = testable.solution - @test integral(testable.integrand, testable.geometry, rule) ≈ sol rtol=rtol + @test integral(testable.integrand, testable.geometry, rule)≈sol rtol=rtol # Callable integrand f = Callable(testable.integrand) - @test integral(f, testable.geometry, rule) ≈ sol rtol=rtol + @test integral(f, testable.geometry, rule)≈sol rtol=rtol # Vector integrand fv(p) = fill(testable.integrand(p), 3) sol_v = fill(testable.solution, 3) - @test integral(fv, testable.geometry, rule) ≈ sol_v rtol=rtol + @test integral(fv, testable.geometry, rule)≈sol_v rtol=rtol else - @test_throws "not supported" integral(testable.integrand, testable.geometry, rule) + @test_throws "not supported" integral( + testable.integrand, testable.geometry, rule) end # if end # for end # function @@ -150,7 +152,7 @@ end # Package and run tests testable = TestableGeometry(integrand, curve, solution) - runtests(testable, SupportStatus(:line); rtol=0.5e-2) + runtests(testable, SupportStatus(:line); rtol = 0.5e-2) end @testitem "Meshes.Box 1D" setup=[Combinations] begin @@ -167,7 +169,7 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) - runtests(testable, SupportStatus(:line); rtol=1e-6) + runtests(testable, SupportStatus(:line); rtol = 1e-6) end @testitem "Meshes.Box 2D" setup=[Combinations] begin @@ -183,7 +185,7 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) - runtests(testable, SupportStatus(:surface); rtol=1e-6) + runtests(testable, SupportStatus(:surface); rtol = 1e-6) end @testitem "Meshes.Box 3D" setup=[Combinations] begin @@ -200,7 +202,7 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) - runtests(testable, SupportStatus(:volume); rtol=1e-6) + runtests(testable, SupportStatus(:volume); rtol = 1e-6) end @testitem "Meshes.Circle" setup=[Combinations] begin @@ -259,7 +261,7 @@ end # Package and run tests testable = TestableGeometry(integrand, cone, solution) - runtests(testable, SupportStatus(:surface); rtol=1e-6) + runtests(testable, SupportStatus(:surface); rtol = 1e-6) end @testitem "Meshes.Cylinder" setup=[Combinations] begin @@ -326,7 +328,7 @@ end # Package and run tests # Tolerances are higher due to `measure` being only an approximation testable = TestableGeometry(integrand, ellipsoid, solution) - runtests(testable, SupportStatus(:surface); rtol=1e-2) + runtests(testable, SupportStatus(:surface); rtol = 1e-2) end @testitem "Meshes.FrustumSurface" setup=[Combinations] begin @@ -358,7 +360,7 @@ end # Package and run tests testable = TestableGeometry(integrand, frustum, solution) - runtests(testable, SupportStatus(:surface); rtol=1e-6) + runtests(testable, SupportStatus(:surface); rtol = 1e-6) end @testitem "Meshes.Hexahedron" setup=[Combinations] begin diff --git a/test/utils.jl b/test/utils.jl index 448d38de..dcddbf3d 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -31,7 +31,7 @@ end @test FiniteDifference().ε ≈ 1e-6 # Test jacobian with wrong number of parametric coordinates - box = Box(Point(0, 0), Point(1, 1)) + box = Box(Point(0, 0), Point(1, 1)) @test_throws ArgumentError jacobian(box, zeros(3)) end From 0d7ab7372e10041d235f842a7070a68b7f6416c1 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 6 Dec 2024 10:38:43 -0500 Subject: [PATCH 37/47] Add framework for planned AutoEnzyme testing --- test/combinations.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index 5ed5bc95..74e0045a 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -89,8 +89,18 @@ else @test_throws "not supported" integral( testable.integrand, testable.geometry, rule) - end # if + end end # for + + #= + iter_diff_methods = ( + (supports.autoenzyme, AutoEnzyme()), + ) + + for (supported, diff_method) in iter_diff_methods + @test integral(testable.integrand, testable.geometry; diff_method=diff_method)≈sol rtol=rtol + end + =# end # function end #testsnippet From 4a512accc57f33feb307a842ff9693f8d37b75d5 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Fri, 6 Dec 2024 16:37:08 -0500 Subject: [PATCH 38/47] Improve code clarity --- test/combinations.jl | 34 +++++++++++++++++++++++----------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 4e015114..23d0c20f 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -1,6 +1,13 @@ -# This section tests: -# - All supported combinations of integral(f, ::Geometry, ::IntegrationAlgorithm) -# - Invalid applications of integral aliases produce a descriptive error +""" +This file includes tests for: +- Supported combinations {::Geometry, ::IntegrationRule} for `integral(f, geometry, rule)`. +- Integrand functions returning a `Unitful.Quantity` or a `Vector{Unitful.Quantity}`. +- Callable objects in place of an integrand function. +- Unsupported combinations produce a descriptive and useful error message. +- The appropriate alias function, e.g. `lineintegral` for 1D geometries, works. +- Invalid applications of integral aliases produce a descriptive and useful error message. +- (Planned) Usage of non-default DifferentiationMethods. +""" #=============================================================================== Test Generation Infrastructure @@ -12,10 +19,11 @@ using MeshIntegrals using Unitful + # Used for testing callable objects as integrand functions struct Callable{F <: Function} f::F end - (c::Callable)(p) = c.f(p) + (c::Callable)(p::Meshes.Point) = c.f(p) # Stores a testable combination struct TestableGeometry{F <: Function, G <: Geometry, U <: Unitful.Quantity} @@ -24,14 +32,18 @@ solution::U end - # Indicates which functions/rules are supported for a particular geometry + # Used to indicate which features are supported for a particular geometry struct SupportStatus + # Alias Functions lineintegral::Bool surfaceintegral::Bool volumeintegral::Bool + # IntegrationRules gausskronrod::Bool gausslegendre::Bool hadaptivecubature::Bool + # DifferentiationMethods + # autoenzyme::Bool end # Shortcut constructor for geometries with typical support structure @@ -56,8 +68,8 @@ function runtests(testable::TestableGeometry, supports::SupportStatus; rtol=sqrt(eps())) # Test alias functions for alias in (lineintegral, surfaceintegral, volumeintegral) - alias_symbol = first(methods(alias)).name - if getfield(supports, alias_symbol) + # if supports.alias + if getfield(supports, first(methods(alias)).name) @test alias(testable.integrand, testable.geometry) ≈ testable.solution rtol=rtol else @test_throws "not supported" alias(testable.integrand, testable.geometry) @@ -177,9 +189,9 @@ end # Integrand & Solution function integrand(p::Meshes.Point) x, y = ustrip.((p.coords.x, p.coords.y)) - (sqrt(a^2 - x^2) + sqrt(a^2 - y^2)) * u"Ω/m^2" + (sqrt(a^2 - x^2) + sqrt(a^2 - y^2)) * u"Ω" end - solution = 2a * (π * a^2 / 4) * u"Ω" + solution = 2a * (π * a^2 / 4) * u"Ω*m^2" # Package and run tests testable = TestableGeometry(integrand, box, solution) @@ -194,9 +206,9 @@ end # Integrand & Solution function integrand(p::Meshes.Point) x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) - (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3" + (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω" end - solution = 3a^2 * (π * a^2 / 4) * u"Ω" + solution = 3a^2 * (π * a^2 / 4) * u"Ω*m^3" # Package and run tests testable = TestableGeometry(integrand, box, solution) From f497e14ce6c290b994e983ebc97f0261decadf0c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 6 Dec 2024 17:56:04 -0500 Subject: [PATCH 39/47] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/combinations.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 7822ebec..2d3a0888 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -66,15 +66,15 @@ This file includes tests for: end function runtests( - testable::TestableGeometry, - supports::SupportStatus; - rtol = sqrt(eps()) + testable::TestableGeometry, + supports::SupportStatus; + rtol = sqrt(eps()) ) # Test alias functions for alias in (lineintegral, surfaceintegral, volumeintegral) # if supports.alias if getfield(supports, first(methods(alias)).name) - @test alias(testable.integrand, testable.geometry) ≈ testable.solution rtol=rtol + @test alias(testable.integrand, testable.geometry)≈testable.solution rtol=rtol else @test_throws "not supported" alias(testable.integrand, testable.geometry) end From ba80123be04fe44f4f8afc841925854282cad973 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 7 Dec 2024 08:56:19 -0500 Subject: [PATCH 40/47] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- test/combinations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 2d3a0888..c8985ec3 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -73,7 +73,7 @@ This file includes tests for: # Test alias functions for alias in (lineintegral, surfaceintegral, volumeintegral) # if supports.alias - if getfield(supports, first(methods(alias)).name) + if getfield(supports, nameof(alias)) @test alias(testable.integrand, testable.geometry)≈testable.solution rtol=rtol else @test_throws "not supported" alias(testable.integrand, testable.geometry) @@ -458,7 +458,7 @@ end # Package and run tests testable_cart = TestableGeometry(integrand, curve_cart, solution) runtests(testable_cart, SupportStatus(:line)) - testable_polar = TestableGeometry(integrand, curve_cart, solution) + testable_polar = TestableGeometry(integrand, curve_polar, solution) runtests(testable_polar, SupportStatus(:line)) end end From 4099736df8bcd65009ea3a01e4abcafa8366ac73 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 7 Dec 2024 21:26:11 -0500 Subject: [PATCH 41/47] Switch from symbols to paramdim [skip ci] Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- test/combinations.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index c8985ec3..44404693 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -47,21 +47,23 @@ This file includes tests for: end # Shortcut constructor for geometries with typical support structure - function SupportStatus(sym::Symbol) - if sym == :line + function SupportStatus(geometry::Geometry) + if paramdim(geometry) == 1 aliases = Bool.((1, 0, 0)) rules = Bool.((1, 1, 1)) return SupportStatus(aliases..., rules...) - elseif sym == :surface + elseif paramdim(geometry) == 2 aliases = Bool.((0, 1, 0)) rules = Bool.((1, 1, 1)) return SupportStatus(aliases..., rules...) - elseif sym == :volume + elseif paramdim(geometry) == 3 aliases = Bool.((0, 0, 1)) rules = Bool.((0, 1, 1)) return SupportStatus(aliases..., rules...) else - error("Unrecognized SupportStatus shortcut $(string(sym))") + aliases = Bool.((0, 0, 0)) + rules = Bool.((0, 1, 1)) + return SupportStatus(aliases..., rules...) end end From e04b75bd7ffecd30948c353555e9716f027fd7bf Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Sat, 7 Dec 2024 21:30:02 -0500 Subject: [PATCH 42/47] Use an optional arg for SupportStatus in runtests --- test/combinations.jl | 48 ++++++++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 44404693..a3e21fa8 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -69,7 +69,7 @@ This file includes tests for: function runtests( testable::TestableGeometry, - supports::SupportStatus; + supports::SupportStatus = SupportStatus(testable.geometry); rtol = sqrt(eps()) ) # Test alias functions @@ -141,7 +141,7 @@ end #testsnippet # Package and run tests testable = TestableGeometry(integrand, ball, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Ball 3D" setup=[Combinations] begin @@ -164,7 +164,7 @@ end # Package and run tests testable = TestableGeometry(integrand, ball, solution) - runtests(testable, SupportStatus(:volume)) + runtests(testable) end @testitem "Meshes.BezierCurve" setup=[Combinations] begin @@ -251,7 +251,7 @@ end # Package and run tests testable = TestableGeometry(integrand, circle, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.Cone" setup=[Combinations] begin @@ -270,7 +270,7 @@ end # Package and run tests testable = TestableGeometry(integrand, cone, solution) - runtests(testable, SupportStatus(:volume)) + runtests(testable) end @testitem "Meshes.ConeSurface" setup=[Combinations] begin @@ -304,7 +304,7 @@ end # Package and run tests testable = TestableGeometry(integrand, cyl, solution) - runtests(testable, SupportStatus(:volume)) + runtests(testable) end @testitem "Meshes.CylinderSurface" setup=[Combinations] begin @@ -319,7 +319,7 @@ end # Package and run tests testable = TestableGeometry(integrand, cyl, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Disk" setup=[Combinations] begin @@ -340,7 +340,7 @@ end # Package and run tests testable = TestableGeometry(integrand, disk, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Ellipsoid" setup=[Combinations] begin @@ -402,7 +402,7 @@ end # Package and run tests testable = TestableGeometry(integrand, hexahedron, solution) - runtests(testable, SupportStatus(:volume)) + runtests(testable) end @testitem "Meshes.Line" setup=[Combinations] begin @@ -420,7 +420,7 @@ end # Package and run tests testable = TestableGeometry(integrand, line, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.ParaboloidSurface" setup=[Combinations] begin @@ -433,7 +433,7 @@ end # Package and run tests testable = TestableGeometry(integrand, parab, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "ParametrizedCurve" setup=[Combinations] begin @@ -459,9 +459,9 @@ end # Package and run tests testable_cart = TestableGeometry(integrand, curve_cart, solution) - runtests(testable_cart, SupportStatus(:line)) + runtests(testable_cart) testable_polar = TestableGeometry(integrand, curve_polar, solution) - runtests(testable_polar, SupportStatus(:line)) + runtests(testable_polar) end end @@ -480,7 +480,7 @@ end # Package and run tests testable = TestableGeometry(integrand, plane, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Quadrangle" setup=[Combinations] begin @@ -498,7 +498,7 @@ end # Package and run tests testable = TestableGeometry(integrand, quadrangle, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Ray" setup=[Combinations] begin @@ -516,7 +516,7 @@ end # Package and run tests testable = TestableGeometry(integrand, ray, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.Ring" setup=[Combinations] begin @@ -536,7 +536,7 @@ end # Package and run tests testable = TestableGeometry(integrand, ring, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.Rope" setup=[Combinations] begin @@ -556,7 +556,7 @@ end # Package and run tests testable = TestableGeometry(integrand, rope, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.Segment" setup=[Combinations] begin @@ -576,7 +576,7 @@ end # Package and run tests testable = TestableGeometry(integrand, segment, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.Sphere 2D" setup=[Combinations] begin @@ -594,7 +594,7 @@ end # Package and run tests testable = TestableGeometry(integrand, sphere, solution) - runtests(testable, SupportStatus(:line)) + runtests(testable) end @testitem "Meshes.Sphere 3D" setup=[Combinations] begin @@ -617,7 +617,7 @@ end # Package and run tests testable = TestableGeometry(integrand, sphere, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Tetrahedron" setup=[Combinations] begin @@ -634,7 +634,7 @@ end # Package and run tests testable = TestableGeometry(integrand, tetrahedron, solution) - runtests(testable, SupportStatus(:volume)) + runtests(testable) end @testitem "Meshes.Torus" setup=[Combinations] begin @@ -649,7 +649,7 @@ end # Package and run tests testable = TestableGeometry(integrand, torus, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end @testitem "Meshes.Triangle" setup=[Combinations] begin @@ -665,5 +665,5 @@ end # Package and run tests testable = TestableGeometry(integrand, triangle, solution) - runtests(testable, SupportStatus(:surface)) + runtests(testable) end From 49473ad54b2cda2d3045281836864e9d5054e9b4 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Sat, 7 Dec 2024 21:38:40 -0500 Subject: [PATCH 43/47] One testsnippet per file --- test/floating_point_types.jl | 11 ++++++++--- test/runtests.jl | 6 ------ test/utils.jl | 16 ++++++++-------- 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/test/floating_point_types.jl b/test/floating_point_types.jl index 57d5a00b..7af87bd5 100644 --- a/test/floating_point_types.jl +++ b/test/floating_point_types.jl @@ -2,14 +2,19 @@ # have expected level of accuracy and are produce results in appropriate type # Base value for atol when integrating with a particular FP type -@testsnippet BaseAtol begin + +@testsnippet FP_Types begin + using LinearAlgebra: norm + using Meshes + using Unitful + baseatol = Dict( Float32 => 0.01f0, BigFloat => BigFloat(0.001) ) end -@testitem "Alternate floating types" setup=[Setup, BaseAtol] begin +@testitem "Alternate floating types" setup=[FP_Types] begin @testset "$FP" for FP in (Float32, BigFloat) # Rectangular volume with unit integrand f = p -> one(FP) @@ -41,7 +46,7 @@ end end end -@testitem "Integral Aliases" setup=[Setup] begin +@testitem "Integral Aliases" setup=[FP_Types] begin f = p -> one(Float32) box1d = Box(Point(fill(0.0f0u"m", 1)...), Point(fill(1.0f0u"m", 1)...)) box2d = Box(Point(fill(0.0f0u"m", 2)...), Point(fill(1.0f0u"m", 2)...)) diff --git a/test/runtests.jl b/test/runtests.jl index 380bff2d..3f981056 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,9 +3,3 @@ using TestItems # For CI, run all tests not marked with the :extended tag @run_package_tests filter=ti -> !(:extended in ti.tags) verbose=true - -@testsnippet Setup begin - using LinearAlgebra: norm - using Meshes - using Unitful -end diff --git a/test/utils.jl b/test/utils.jl index dcddbf3d..947688a7 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,7 +1,11 @@ -@testitem "Utilities" setup=[Setup] begin +@testsnippet Utils begin using LinearAlgebra: norm - using MeshIntegrals: _units, _zeros, _ones + using Meshes + using MeshIntegrals: _default_diff_method, _parametric, _units, _zeros, _ones + using Unitful +end +@testitem "Utilities" setup=[Utils] begin # _KVector v = Meshes.Vec(3, 4) @test norm(MeshIntegrals._KVector(v)) ≈ 5.0u"m" @@ -19,9 +23,7 @@ @test _ones(Float32, 2) == (1.0f0, 1.0f0) end -@testitem "Differentiation" setup=[Setup] begin - using MeshIntegrals: _default_diff_method - +@testitem "Differentiation" setup=[Utils] begin # _default_diff_method sphere = Sphere(Point(0, 0, 0), 1.0) @test _default_diff_method(Meshes.Sphere) isa FiniteDifference @@ -35,9 +37,7 @@ end @test_throws ArgumentError jacobian(box, zeros(3)) end -@testitem "_ParametricGeometry" setup=[Setup] begin - using MeshIntegrals: _parametric - +@testitem "_ParametricGeometry" setup=[Utils] begin pt_n = Point(0, 3, 0) pt_w = Point(-7, 0, 0) pt_e = Point(8, 0, 0) From 3f8b9f0ffc7e60135fc8b2910c83a6efc04a42c1 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Sat, 7 Dec 2024 21:42:44 -0500 Subject: [PATCH 44/47] Update some missed constructor changes --- test/combinations.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index a3e21fa8..7190d65d 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -180,7 +180,7 @@ end # Package and run tests testable = TestableGeometry(integrand, curve, solution) - runtests(testable, SupportStatus(:line); rtol = 0.5e-2) + runtests(testable; rtol = 0.5e-2) end @testitem "Meshes.Box 1D" setup=[Combinations] begin @@ -197,7 +197,7 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) - runtests(testable, SupportStatus(:line); rtol = 1e-6) + runtests(testable; rtol = 1e-6) end @testitem "Meshes.Box 2D" setup=[Combinations] begin @@ -213,7 +213,7 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) - runtests(testable, SupportStatus(:surface); rtol = 1e-6) + runtests(testable; rtol = 1e-6) end @testitem "Meshes.Box 3D" setup=[Combinations] begin @@ -230,7 +230,7 @@ end # Package and run tests testable = TestableGeometry(integrand, box, solution) - runtests(testable, SupportStatus(:volume); rtol = 1e-6) + runtests(testable; rtol = 1e-6) end @testitem "Meshes.Circle" setup=[Combinations] begin @@ -289,7 +289,7 @@ end # Package and run tests testable = TestableGeometry(integrand, cone, solution) - runtests(testable, SupportStatus(:surface); rtol = 1e-6) + runtests(testable; rtol = 1e-6) end @testitem "Meshes.Cylinder" setup=[Combinations] begin @@ -356,7 +356,7 @@ end # Package and run tests # Tolerances are higher due to `measure` being only an approximation testable = TestableGeometry(integrand, ellipsoid, solution) - runtests(testable, SupportStatus(:surface); rtol = 1e-2) + runtests(testable; rtol = 1e-2) end @testitem "Meshes.FrustumSurface" setup=[Combinations] begin @@ -388,7 +388,7 @@ end # Package and run tests testable = TestableGeometry(integrand, frustum, solution) - runtests(testable, SupportStatus(:surface); rtol = 1e-6) + runtests(testable; rtol = 1e-6) end @testitem "Meshes.Hexahedron" setup=[Combinations] begin From 482caf46cbaad542276d76329f0ffd0e687cbc38 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Sat, 7 Dec 2024 21:57:01 -0500 Subject: [PATCH 45/47] Restore Box 4D tests --- test/combinations.jl | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index 7190d65d..fa0e724e 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -190,10 +190,10 @@ end # Integrand & Solution function integrand(p::Meshes.Point) - t = ustrip(p.coords.x) - sqrt(a^2 - t^2) * u"Ω" + x₁ = ustrip.((to(p))) + √(a^2 - x₁^2) * u"A" end - solution = π * a^2 / 4 * u"Ω*m" + solution = π * a^2 / 4 * u"A*m" # Package and run tests testable = TestableGeometry(integrand, box, solution) @@ -206,10 +206,10 @@ end # Integrand & Solution function integrand(p::Meshes.Point) - x, y = ustrip.((p.coords.x, p.coords.y)) - (sqrt(a^2 - x^2) + sqrt(a^2 - y^2)) * u"Ω" + x₁, x₂ = ustrip.((to(p))) + (√(a^2 - x₁^2) + √(a^2 - x₂^2)) * u"A" end - solution = 2a * (π * a^2 / 4) * u"Ω*m^2" + solution = 2a * (π * a^2 / 4) * u"A*m^2" # Package and run tests testable = TestableGeometry(integrand, box, solution) @@ -223,10 +223,27 @@ end # Integrand & Solution function integrand(p::Meshes.Point) - x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) - (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω" + x₁, x₂, x₃ = ustrip.((to(p))) + (√(a^2 - x₁^2) + √(a^2 - x₂^2) + √(a^2 - x₃^2)) * u"A" + end + solution = 3a^2 * (π * a^2 / 4) * u"A*m^3" + + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable; rtol = 1e-6) +end + +@testitem "Meshes.Box 4D" setup=[Combinations] begin + # Geometry + a = π + box = Box(Point(0, 0, 0, 0), Point(a, a, a, a)) + + # Integrand & Solution + function integrand(p::Meshes.Point) + x₁, x₂, x₃, x₄ = ustrip.((to(p))) + (√(a^2 - x₁^2) + √(a^2 - x₂^2) + √(a^2 - x₃^2) + √(a^2 - x₄^2)) * u"A" end - solution = 3a^2 * (π * a^2 / 4) * u"Ω*m^3" + solution = 4a^3 * (π * a^2 / 4) * u"A*m^4" # Package and run tests testable = TestableGeometry(integrand, box, solution) From dcbfae471b0c31154d17e5e16953ea97e12c0105 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Sat, 7 Dec 2024 22:01:02 -0500 Subject: [PATCH 46/47] Restore extended tag --- test/combinations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index fa0e724e..d7757e08 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -233,7 +233,7 @@ end runtests(testable; rtol = 1e-6) end -@testitem "Meshes.Box 4D" setup=[Combinations] begin +@testitem "Meshes.Box 4D" tags=[:extended] setup=[Combinations] begin # Geometry a = π box = Box(Point(0, 0, 0, 0), Point(a, a, a, a)) From 9c05943dee1c4bf32ec5996d061ccbf6d731d7fd Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 7 Dec 2024 22:21:25 -0500 Subject: [PATCH 47/47] Reduce dims to single coordinate --- test/combinations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/combinations.jl b/test/combinations.jl index d7757e08..b53420bd 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -190,7 +190,7 @@ end # Integrand & Solution function integrand(p::Meshes.Point) - x₁ = ustrip.((to(p))) + x₁ = only(ustrip.((to(p)))) √(a^2 - x₁^2) * u"A" end solution = π * a^2 / 4 * u"A*m"