From c7c0a47fc475b96599c8ed8337f313f4de732830 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 8 Dec 2024 07:18:30 -0500 Subject: [PATCH] Standardization of combinations tests (#151) * Initial attempt at new test generation * Generate alias tests in a loop * Loop over all integration rules * Implement for Circle * Code cleanup * Bugfix * Implement for Cone * Implement for Cylinder * Implement for CylinderSurface and Disk * Formatting * Implement for Line, ParaboloidSurface, and ParametrizedCurve * Update spelling * Update setup * Implement for Plane, Quadrangle, and Ray * Update integrand name * Implement for Triangle * Implement for Torus * Implement Spheres * Fix typos * Implement for Rope and Segment * Fix units * Implement for Ring * Absorb integrand parameters * Code re-org and remove Box-4D * Code cleanup and add Unitful integrands * Code re-org in Hexahedron * Typo fix * Bugfixes * Add rtol kwarg to runtests * Implement for BezierCurve and Box-1D * Implement for Box 2D and 3D * Implement ConeSurface and Ellipsoid * Implement for FrustumSurface and Hexahedron * Move jacobian tests * Bugfix * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Add framework for planned AutoEnzyme testing * Improve code clarity * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Switch from symbols to paramdim [skip ci] Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Use an optional arg for SupportStatus in runtests * One testsnippet per file * Update some missed constructor changes * Restore Box 4D tests * Restore extended tag * Reduce dims to single coordinate --------- Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- test/combinations.jl | 1108 ++++++++++++++-------------------- test/floating_point_types.jl | 11 +- test/runtests.jl | 6 - test/utils.jl | 20 +- 4 files changed, 461 insertions(+), 684 deletions(-) diff --git a/test/combinations.jl b/test/combinations.jl index a6d63d3d..b53420bd 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -1,429 +1,384 @@ -# 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 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 +===============================================================================# + +@testsnippet Combinations begin + using LinearAlgebra: norm + using Meshes + using MeshIntegrals + using Unitful + + # Used for testing callable objects as integrand functions + struct Callable{F <: Function} + f::F + end + (c::Callable)(p::Meshes.Point) = c.f(p) + + # Stores a testable combination + struct TestableGeometry{F <: Function, G <: Geometry, U <: Unitful.Quantity} + integrand::F + geometry::G + solution::U + end + + # 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 + function SupportStatus(geometry::Geometry) + if paramdim(geometry) == 1 + aliases = Bool.((1, 0, 0)) + rules = Bool.((1, 1, 1)) + return SupportStatus(aliases..., rules...) + elseif paramdim(geometry) == 2 + aliases = Bool.((0, 1, 0)) + rules = Bool.((1, 1, 1)) + return SupportStatus(aliases..., rules...) + elseif paramdim(geometry) == 3 + aliases = Bool.((0, 0, 1)) + rules = Bool.((0, 1, 1)) + return SupportStatus(aliases..., rules...) + else + aliases = Bool.((0, 0, 0)) + rules = Bool.((0, 1, 1)) + return SupportStatus(aliases..., rules...) + end + end + + function runtests( + testable::TestableGeometry, + supports::SupportStatus = SupportStatus(testable.geometry); + rtol = sqrt(eps()) + ) + # Test alias functions + for alias in (lineintegral, surfaceintegral, volumeintegral) + # if supports.alias + 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) + end + end -@testitem "Meshes.Ball 2D" setup=[Setup] begin + iter_rules = ( + (supports.gausskronrod, GaussKronrod()), + (supports.gausslegendre, GaussLegendre(100)), + (supports.hadaptivecubature, HAdaptiveCubature()) + ) + + # Test rules + for (supported, rule) in iter_rules + if supported + # Scalar integrand + sol = testable.solution + @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 + + # 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 + else + f = testable.integrand + geometry = testable.geometry + @test_throws "not supported" integral(f, geometry, rule) + 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 + +#=============================================================================== + Create and Test Geometries +===============================================================================# + +@testitem "Meshes.Ball 2D" setup=[Combinations] begin + # Geometry origin = Point(0, 0) radius = 2.8 ball = Ball(origin, radius) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(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 = (π - π * 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 - - # 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) + solution = (π - π * exp(-radius^2)) * u"A*m^2" + + # Package and run tests + testable = TestableGeometry(integrand, ball, solution) + runtests(testable) end -@testitem "Meshes.Ball 3D" setup=[Setup] begin +@testitem "Meshes.Ball 3D" setup=[Combinations] begin using SpecialFunctions: erf + # Geometry center = Point(1, 2, 3) radius = 2.8u"m" + r = ustrip(u"m", radius) ball = Ball(center, radius) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(p::Meshes.Point) offset = p - center ur = norm(offset) r = ustrip(u"m", ur) - exp(-r^2) + exp(-r^2) * u"A" end - fv(p) = fill(f(p), 3) + solution = (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"A*m^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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, ball) - @test_throws "not supported" surfaceintegral(f, ball) - @test volumeintegral(f, ball) ≈ sol + # Package and run tests + testable = TestableGeometry(integrand, ball, solution) + runtests(testable) end -@testitem "Meshes.BezierCurve" setup=[Setup] begin - curve = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)]) +@testitem "Meshes.BezierCurve" setup=[Combinations] begin + # Geometry + curve = BezierCurve([Point(t, sin(t), 0) for t in range(-π, π, length = 361)]) - function f(p::P) where {P <: Meshes.Point} + # Integrand + function integrand(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) - - # 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 - - # 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 - - # 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) - - # Check Bezier-specific jacobian bounds - @test_throws DomainError jacobian(curve, [1.1]) + solution = 2π * u"Ω*m" + + # Package and run tests + testable = TestableGeometry(integrand, curve, solution) + runtests(testable; rtol = 0.5e-2) end -@testitem "Meshes.Box 1D" setup=[Setup] begin +@testitem "Meshes.Box 1D" setup=[Combinations] begin + # Geometry a = π box = Box(Point(0), Point(a)) - struct Integrand end - fc = Integrand() - - function (::Integrand)(p::Meshes.Point) - t = ustrip(p.coords.x) - sqrt(a^2 - t^2) * u"Ω/m" + # Integrand & Solution + function integrand(p::Meshes.Point) + x₁ = only(ustrip.((to(p)))) + √(a^2 - x₁^2) * u"A" end - fcv(p) = fill(fc(p), 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 - - # 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 - - # Integral aliases - @test lineintegral(fc, box) ≈ sol - @test_throws "not supported" surfaceintegral(fc, box) - @test_throws "not supported" volumeintegral(fc, box) + solution = π * a^2 / 4 * u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable; 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::P) where {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" + # Integrand & Solution + function integrand(p::Meshes.Point) + x₁, x₂ = ustrip.((to(p))) + (√(a^2 - x₁^2) + √(a^2 - x₂^2)) * u"A" end - fv(p) = fill(f(p), 3) - - # 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) - - # Test jacobian with wrong number of parametric coordinates - @test_throws ArgumentError jacobian(box, zeros(3)) + solution = 2a * (π * a^2 / 4) * u"A*m^2" + + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable; rtol = 1e-6) 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)) - function f(p::P) where {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" + # Integrand & Solution + function integrand(p::Meshes.Point) + x₁, x₂, x₃ = ustrip.((to(p))) + (√(a^2 - x₁^2) + √(a^2 - x₂^2) + √(a^2 - x₃^2)) * u"A" end - fv(p) = fill(f(p), 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test_throws "not supported" surfaceintegral(f, box) - @test volumeintegral(f, box) ≈ sol + 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" tags=[:extended] setup=[Setup] begin +@testitem "Meshes.Box 4D" tags=[:extended] setup=[Combinations] begin + # Geometry a = π box = Box(Point(0, 0, 0, 0), Point(a, a, a, a)) - function f(p::P) where {P <: Meshes.Point} - x1, x2, x3, x4 = ustrip.(to(p).coords) - σ(x) = sqrt(a^2 - x^2) - (σ(x1) + σ(x2) + σ(x3) + σ(x4)) * u"Ω/m^4" + # 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 - 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) + solution = 4a^3 * (π * a^2 / 4) * u"A*m^4" + + # Package and run tests + testable = TestableGeometry(integrand, box, solution) + runtests(testable; rtol = 1e-6) 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 & Solution + function integrand(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) - - # 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 - - # 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) + solution = 2π * radius * exp(-radius^2) * u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, circle, solution) + runtests(testable) 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 + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = (π * r^2 * h / 3) * u"A" - # 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 - - # 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 - - # 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) 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) apex = Point(0.0u"m", 0.0u"m", h) cone = ConeSurface(base, apex) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - - _area_cone_rightcircular(h, r) = (π * r^2) + (π * r * hypot(h, r)) + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = ((π * r^2) + (π * r * hypot(h, r))) * u"A" - # 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 - - # 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; rtol = 1e-6) 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 - - # 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, cyl) - @test_throws "not supported" surfaceintegral(f, cyl) - @test volumeintegral(f, cyl) ≈ sol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(cyl) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, cyl, solution) + runtests(testable) 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, cyl) - @test surfaceintegral(f, cyl) ≈ sol - @test_throws "not supported" volumeintegral(f, cyl) + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(cyl) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, cyl, solution) + runtests(testable) 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::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) - - # 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) + solution = (π - π * exp(-radius^2)) * u"A*m^2" + + # Package and run tests + testable = TestableGeometry(integrand, disk, solution) + runtests(testable) 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) - f(p) = 1.0 - fv(p) = fill(f(p), 3) + # Integrand & Solution + 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 - 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 - - # 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; 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" @@ -436,477 +391,296 @@ end disk_top = Disk(plane_top, r_top) frustum = FrustumSurface(disk_bot, disk_top) - f(p) = 1.0 - fv(p) = fill(f(p), 3) - + # Integrand & Solution + integrand(p) = 1.0u"A" _area_base(r) = π * r^2 _area_cone_walls(h, r) = π * r * hypot(h, r) - - # Scalar integrand - 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_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 - @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 + + # Package and run tests + testable = TestableGeometry(integrand, frustum, solution) + runtests(testable; 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)) - f(p) = 1.0 - fv(p) = fill(f(p), 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, hexahedron) - @test_throws "not supported" surfaceintegral(f, hexahedron) - @test volumeintegral(f, hexahedron) ≈ sol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(hexahedron) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, hexahedron, solution) + runtests(testable) 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") +@testitem "Meshes.Line" setup=[Combinations] begin + # Geometry + a = Point(0, 0, 0) + b = Point(1, 1, 1) line = Line(a, b) - function f(p::P) where {P <: Meshes.Point} + # Integrand & solution + function integrand(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 - - # 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) + solution = sqrt(π) * u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, line, solution) + runtests(testable) 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, parab) - @test surfaceintegral(f, parab) ≈ sol - @test_throws "not supported" volumeintegral(f, parab) + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(parab) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, parab, solution) + runtests(testable) 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 # 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π)) curve_polar = ParametrizedCurve(t -> Point(Polar(radius, t)), (0.0, 2π)) - function f(p::P) where {P <: Meshes.Point} + # Integrand & Solution + function integrand(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) + testable_polar = TestableGeometry(integrand, curve_polar, solution) + runtests(testable_polar) 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::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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, plane) - @test surfaceintegral(f, plane) ≈ sol - @test_throws "not supported" volumeintegral(f, plane) + solution = π * u"A*m^2" + + # Package and run tests + testable = TestableGeometry(integrand, plane, solution) + runtests(testable) 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::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 = 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) + solution = 0.5 * π * erf(1)^2 * u"A*m^2" + + # Package and run tests + testable = TestableGeometry(integrand, quadrangle, solution) + runtests(testable) 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::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(π) / 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) + solution = sqrt(π) / 2 * u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, ray, solution) + runtests(testable) 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::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 - - # 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) + solution = 14.0u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, ring, solution) + runtests(testable) 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") - rope = Rope(pt_a, pt_b, pt_c, pt_d) +@testitem "Meshes.Rope" setup=[Combinations] begin + # Geometry + 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) - 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::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 - - # Integral aliases - @test lineintegral(f, rope) ≈ sol - @test_throws "not supported" surfaceintegral(f, rope) - @test_throws "not supported" volumeintegral(f, rope) + solution = 7.0u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, rope, solution) + runtests(testable) 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 - - # Integral aliases - @test lineintegral(f, segment) ≈ sol - @test_throws "not supported" surfaceintegral(f, segment) - @test_throws "not supported" volumeintegral(f, segment) + solution = ((a - b) / (log(a) - log(b))) * u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, segment, solution) + runtests(testable) 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::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 = 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 - - # Integral aliases - @test lineintegral(f, sphere) ≈ sol - @test_throws "not supported" surfaceintegral(f, sphere) - @test_throws "not supported" volumeintegral(f, sphere) + solution = 2π * radius * exp(-radius^2) * u"A*m" + + # Package and run tests + testable = TestableGeometry(integrand, sphere, solution) + runtests(testable) 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::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) - - # 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) + solution = ((2π * radius^2) + (4π / 3 * radius^2)) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, sphere, solution) + runtests(testable) 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, tetrahedron) - @test_throws "not supported" surfaceintegral(f, tetrahedron) - @test volumeintegral(f, tetrahedron, HAdaptiveCubature()) ≈ sol + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(tetrahedron) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, tetrahedron, solution) + runtests(testable) 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 - - # Integral aliases - @test_throws "not supported" lineintegral(f, torus) - @test surfaceintegral(f, torus) ≈ sol - @test_throws "not supported" volumeintegral(f, torus) + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(torus) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, torus, solution) + runtests(testable) 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 - 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) + # Integrand & Solution + integrand(p) = 1.0u"A" + solution = Meshes.measure(triangle) * u"A" + + # Package and run tests + testable = TestableGeometry(integrand, triangle, solution) + runtests(testable) end 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 50cb4ea1..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 "DifferentiationMethod" 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 @@ -29,11 +31,13 @@ end # FiniteDifference @test FiniteDifference().ε ≈ 1e-6 -end -@testitem "_ParametricGeometry" setup=[Setup] begin - using MeshIntegrals: _parametric + # Test jacobian with wrong number of parametric coordinates + box = Box(Point(0, 0), Point(1, 1)) + @test_throws ArgumentError jacobian(box, zeros(3)) +end +@testitem "_ParametricGeometry" setup=[Utils] begin pt_n = Point(0, 3, 0) pt_w = Point(-7, 0, 0) pt_e = Point(8, 0, 0)