Skip to content

Commit

Permalink
Merge f45d7bb into 50268f4
Browse files Browse the repository at this point in the history
  • Loading branch information
@mikeingold
mikeingold authored Sep 8, 2024
2 parents 50268f4 + f45d7bb commit 8b0bcbe
Show file tree
Hide file tree
Showing 2 changed files with 159 additions and 115 deletions.
4 changes: 2 additions & 2 deletions src/integral_aliases.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ function surfaceintegral(
Dim = Meshes.paramdim(geometry)

if Dim == 2
return integral(f, geometry, GaussKronrod())
return integral(f, geometry, HAdaptiveCubature())
else
error("Performing a surface integral on a geometry with $Dim parametric dimensions not supported.")
end
Expand Down Expand Up @@ -145,7 +145,7 @@ function volumeintegral(
Dim = Meshes.paramdim(geometry)

if Dim == 3
return integral(f, geometry, GaussKronrod())
return integral(f, geometry, HAdaptiveCubature())
else
error("Performing a volume integral on a geometry with $Dim parametric dimensions not supported.")
end
Expand Down
270 changes: 157 additions & 113 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,22 @@
using Aqua
using MeshIntegrals
using Meshes
using MeshIntegrals
using Test
using Unitful


################################################################################
# Aqua.jl Tests
################################################################################

@testset "Aqua.jl" begin
# As of v0.11.4:
# - Ambiguities check disabled since it fails due to upstream findings
# - Verified that no ambiguities exist within MeshIntegrals.jl
Aqua.test_all(MeshIntegrals; ambiguities=false)
end


################################################################################
# Infrastructure
################################################################################
Expand Down Expand Up @@ -142,127 +155,77 @@ end
@testset "Float64 Geometries" begin
map(autotest, SUPPORT_MATRIX(Float64))
end
end

@testset "Float32 Geometries" begin
# TODO temp disabled, see Issue #33
#map(autotest, SUPPORT_MATRIX(Float64))
end

# Custom tests for Line (no measure available for reference)
@testset "Meshes.Line" begin
line = Line(pt_e(Float64), pt_w(Float64))

function f(p::P) where {P<:Meshes.Point}
x = ustrip(u"m", p.coords.x)
y = ustrip(u"m", p.coords.y)
z = ustrip(u"m", p.coords.z)
exp(-x^2)
end
fv(p) = fill(f(p),3)

@test integral(f, line, GaussLegendre(100)) ≈ sqrt(π)*u"m"
@test integral(f, line, GaussKronrod()) ≈ sqrt(π)*u"m"
@test integral(f, line, HAdaptiveCubature()) ≈ sqrt(π)*u"m"

@test integral(fv, line, GaussLegendre(100)) ≈ fill(sqrt(π)*u"m",3)
@test integral(fv, line, GaussKronrod()) ≈ fill(sqrt(π)*u"m",3)
@test integral(fv, line, HAdaptiveCubature()) ≈ fill(sqrt(π)*u"m",3)
end

# Custom tests for Ray (no measure available for reference)
@testset "Meshes.Ray" begin
ray = Ray(origin3d(Float64), ẑ(Float64))

function f(p::P) where {P<:Meshes.Point}
x = ustrip(u"m", p.coords.x)
y = ustrip(u"m", p.coords.y)
z = ustrip(u"m", p.coords.z)
2 * exp(-z^2)
end
fv(p) = fill(f(p),3)

@test integral(f, ray, GaussLegendre(100)) ≈ sqrt(π)*u"m"
@test integral(f, ray, GaussKronrod()) ≈ sqrt(π)*u"m"
@test integral(f, ray, HAdaptiveCubature()) ≈ sqrt(π)*u"m"

@test integral(fv, ray, GaussLegendre(100)) ≈ fill(sqrt(π)*u"m",3)
@test integral(fv, ray, GaussKronrod()) ≈ fill(sqrt(π)*u"m",3)
@test integral(fv, ray, HAdaptiveCubature()) ≈ fill(sqrt(π)*u"m",3)
end

# Custom tests for Plane (no measure available for reference)
@testset "Meshes.Plane" begin
plane = Plane(origin3d(Float64), ẑ(Float64))

function f(p::P) where {P<:Meshes.Point}
x = ustrip(u"m", p.coords.x)
y = ustrip(u"m", p.coords.y)
z = ustrip(u"m", p.coords.z)
exp(-x^2 - y^2)
end
fv(p) = fill(f(p),3)

@test integral(f, plane, GaussLegendre(100)) ≈ π*u"m^2"
@test integral(f, plane, GaussKronrod()) ≈ π*u"m^2"
@test integral(f, plane, HAdaptiveCubature()) ≈ π*u"m^2"

################################################################################
# New Tests
################################################################################

@test integral(fv, plane, GaussLegendre(100)) ≈ fill(π*u"m^2",3)
@test integral(fv, plane, GaussKronrod()) ≈ fill(π*u"m^2",3)
@test integral(fv, plane, HAdaptiveCubature()) ≈ fill(π*u"m^2",3)
end
@testset "New Independent Tests" begin

# Custom tests for Cone
@testset "Meshes.Cone" begin
T = Float64

cone_r = T(2.5)
cone_h = T(2.5)

cone = let
base = Disk(plane_xy(T), cone_r)
Cone(base, Point(0, 0, cone_h))
end

f(p) = T(1)
r = 2.5u"m"
h = 2.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) = T(π) * r^2 * h / 3
cone_volume = _volume_cone_rightcircular(cone_r * u"m", cone_h * u"m")
_volume_cone_rightcircular(h, r) = π * r^2 * h / 3

@test integral(f, cone, GaussLegendre(100)) ≈ cone_volume
# 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()) ≈ cone_volume
@test integral(f, cone, HAdaptiveCubature()) ≈ sol

@test integral(fv, cone, GaussLegendre(100)) ≈ fill(cone_volume, 3)
# 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()) ≈ fill(cone_volume, 3)
@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
end

# Custom tests for ConeSurface
@testset "Meshes.ConeSurface" begin
T = Float64

cone_r = T(2.5)
cone_h = T(2.5)

cone = let
base = Disk(plane_xy(T), cone_r)
ConeSurface(base, Point(0, 0, cone_h))
end

f(p) = T(1)
r = 2.5u"m"
h = 2.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) = T(π) * r^2 + T(π) * r * hypot(h, r)
cone_area = _area_cone_rightcircular(cone_r * u"m", cone_h * u"m")
_area_cone_rightcircular(h, r) = (π * r^2) + (π * r * hypot(h, r))

# Scalar integrand
sol = _area_cone_rightcircular(h, r)
@test integral(f, cone, GaussLegendre(100)) ≈ sol
@test integral(f, cone, GaussKronrod()) ≈ sol
@test integral(f, cone, HAdaptiveCubature()) ≈ sol

@test integral(f, cone, GaussLegendre(100)) ≈ cone_area
@test integral(f, cone, GaussKronrod()) ≈ cone_area
@test integral(f, cone, HAdaptiveCubature()) ≈ cone_area
# Vector integrand
vsol = fill(sol, 3)
@test integral(fv, cone, GaussLegendre(100)) ≈ vsol
@test integral(fv, cone, GaussKronrod()) ≈ vsol
@test integral(fv, cone, HAdaptiveCubature()) ≈ vsol

@test integral(fv, cone, GaussLegendre(100)) ≈ fill(cone_area, 3)
@test integral(fv, cone, GaussKronrod()) ≈ fill(cone_area, 3)
@test integral(fv, cone, HAdaptiveCubature()) ≈ fill(cone_area, 3)
# Integral aliases
@test_throws "not supported" lineintegral(f, cone)
@test surfaceintegral(f, cone) ≈ sol
@test_throws "not supported" volumeintegral(f, cone)
end

#= DISABLED FrustumSurface testing due to long run times and seemingly-incorrect results
Expand Down Expand Up @@ -302,14 +265,95 @@ end
@test integral(fv, frustum, HAdaptiveCubature()) ≈ fill(frustum_area, 3)
end
=#
end

################################################################################
# Aqua.jl Tests
################################################################################
@testset "Meshes.Line" begin
a = Point(0.0u"m", 0.0u"m", 0.0u"m")
b = Point(1.0u"m", 1.0u"m", 1.0u"m")
line = Line(a, b)

function f(p::P) where {P<:Meshes.Point}
ur = hypot(p.coords.x, p.coords.y, p.coords.z)
r = ustrip(u"m", ur)
exp(-r^2)
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)
end

@testset "Meshes.Plane" begin
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}
ur = hypot(p.coords.x, p.coords.y, p.coords.z)
r = ustrip(u"m", ur)
exp(-r^2)
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)
end

@testset "Meshes.Ray" begin
a = Point(0.0u"m", 0.0u"m", 0.0u"m")
v = Vec(1.0u"m", 1.0u"m", 1.0u"m")
ray = Ray(a, v)

function f(p::P) where {P<:Meshes.Point}
ur = hypot(p.coords.x, p.coords.y, p.coords.z)
r = ustrip(u"m", ur)
exp(-r^2)
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)
end

@testset "Aqua.jl" begin
# As of v0.11.4 -- Ambiguities check disabled since it fails due to upstream findings
# Verified that no ambiguities exist within MeshIntegrals.jl
Aqua.test_all(MeshIntegrals; ambiguities=false)
end

0 comments on commit 8b0bcbe

Please sign in to comment.