From 1fbd4dde616774f1fa5b3bf438c42760e4a01258 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 25 Feb 2024 11:41:06 -0500 Subject: [PATCH 01/10] Transform line surface and volume functions into dim-checked aliases for integral --- src/integral.jl | 37 +++++++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 4a37cd9e..930087cd 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -56,6 +56,12 @@ end # Integral Function API ################################################################################ +""" + integral(f, geometry, algorithm::IntegrationAlgorithm) + +Numerically integrate a given function `f(::Point)` over the domain defined by +a `geometry` using a particular `integration algorithm`. +""" function integral end """ @@ -67,12 +73,19 @@ using a particular `integration algorithm`. Algorithm types available: - GaussKronrod - GaussLegendre +- HAdaptiveCubature """ function lineintegral( f::F, - geometry::G + geometry::G, + settings::IntegrationAlgorithm=GaussKronrod() ) where {F<:Function, G<:Meshes.Geometry} - return lineintegral(f, geometry, GaussKronrod()) + dim = paramdim(geometry) + if dim == 1 + return integral(f, geometry, settings) + else + error("Unable to perform line integral on a geometry with $dim parametric dimensions.") + end end """ @@ -83,9 +96,15 @@ using a particular `integration algorithm`. """ function surfaceintegral( f::F, - geometry::G + geometry::G, + settings::IntegrationAlgorithm=HAdaptiveCubature() ) where {F<:Function, G<:Meshes.Geometry} - return surfaceintegral(f, geometry, HAdaptiveCubature()) + dim = paramdim(geometry) + if dim == 2 + return integral(f, geometry, settings) + else + error("Unable to perform surface integral on a geometry with $dim parametric dimensions.") + end end """ @@ -97,7 +116,13 @@ Numerically integrate a given function `f(::Point)` throughout a volumetric """ function volumeintegral( f::F, - geometry::G + geometry::G, + settings::IntegrationAlgorithm=HAdaptiveCubature() ) where {F<:Function, G<:Meshes.Geometry} - return volumeintegral(f, geometry, HAdaptiveCubature()) + dim = paramdim(geometry) + if dim == 3 + return integral(f, geometry, settings) + else + error("Unable to perform surface integral on a geometry with $dim parametric dimensions.") + end end From 2273df42d529f86d316beb646ea58ccb74de4451 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 25 Feb 2024 11:43:18 -0500 Subject: [PATCH 02/10] Tweak arg types --- src/integral.jl | 12 ++++----- src/integral_surface.jl | 59 +++++++++++++++++++++++------------------ 2 files changed, 39 insertions(+), 32 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 930087cd..ea0a18dc 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -78,8 +78,8 @@ Algorithm types available: function lineintegral( f::F, geometry::G, - settings::IntegrationAlgorithm=GaussKronrod() -) where {F<:Function, G<:Meshes.Geometry} + settings::I=GaussKronrod() +) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm} dim = paramdim(geometry) if dim == 1 return integral(f, geometry, settings) @@ -97,8 +97,8 @@ using a particular `integration algorithm`. function surfaceintegral( f::F, geometry::G, - settings::IntegrationAlgorithm=HAdaptiveCubature() -) where {F<:Function, G<:Meshes.Geometry} + settings::I=HAdaptiveCubature() +) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm} dim = paramdim(geometry) if dim == 2 return integral(f, geometry, settings) @@ -117,8 +117,8 @@ Numerically integrate a given function `f(::Point)` throughout a volumetric function volumeintegral( f::F, geometry::G, - settings::IntegrationAlgorithm=HAdaptiveCubature() -) where {F<:Function, G<:Meshes.Geometry} + settings::I=HAdaptiveCubature() +) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm} dim = paramdim(geometry) if dim == 3 return integral(f, geometry, settings) diff --git a/src/integral_surface.jl b/src/integral_surface.jl index dcf84e9a..fa034146 100644 --- a/src/integral_surface.jl +++ b/src/integral_surface.jl @@ -1,8 +1,16 @@ +function lineintegral(f, geometry, settings) + dim = paramdim(geometry) + if == 1 + return integral(f, geometry, settings) + else + error("Unable to perform line integral on a geometry with $dim parametric dimensions.") +end + ################################################################################ # Gauss-Legendre ################################################################################ -function surfaceintegral( +function integral( f, disk::Meshes.Ball{2,T}, settings::GaussLegendre @@ -28,7 +36,7 @@ function surfaceintegral( return (0.25 * area(disk)) .* sum(g, zip(wws,xxs)) end -function surfaceintegral( +function integral( f, box::Meshes.Box{2,T}, settings::GaussLegendre @@ -54,7 +62,7 @@ function surfaceintegral( return 0.25 * area(box) .* sum(g, zip(wws,xxs)) end -function surfaceintegral( +function integral( f, disk::Meshes.Disk{T}, settings::GaussLegendre @@ -84,14 +92,14 @@ function surfaceintegral( end """ - surfaceintegral(f, triangle::Meshes.Triangle, ::GaussLegendre) + integral(f, triangle::Meshes.Triangle, ::GaussLegendre) -Like [`surfaceintegral`](@ref) but integrates over the surface of a `triangle` +Like [`integral`](@ref) but integrates over the surface of a `triangle` by transforming the triangle into a polar-barycentric coordinate system and using a Gauss-Legendre quadrature rule along each barycentric dimension of the triangle. """ -function surfaceintegral( +function integral( f, triangle::Meshes.Ngon{3,Dim,T}, settings::GaussLegendre @@ -129,7 +137,7 @@ function surfaceintegral( return (π/4) * area(triangle) .* sum(integrand, zip(wws,xxs)) end -function surfaceintegral( +function integral( f, sphere::Meshes.Sphere{3,T}, settings::GaussLegendre @@ -157,7 +165,7 @@ end # Gauss-Kronrod ################################################################################ -function surfaceintegral( +function integral( f, disk::Meshes.Ball{2,T}, settings::GaussKronrod @@ -177,7 +185,7 @@ function surfaceintegral( return (2π * disk.radius) .* outerintegral end -function surfaceintegral( +function integral( f, box::Meshes.Box{2,T}, settings::GaussKronrod @@ -193,7 +201,7 @@ function surfaceintegral( return area(box) .* outerintegral end -function surfaceintegral( +function integral( f, cyl::Meshes.CylinderSurface{T}, settings::GaussKronrod @@ -228,7 +236,7 @@ function surfaceintegral( return sides + top + bottom end -function surfaceintegral( +function integral( f, disk::Meshes.Disk{T}, settings::GaussKronrod @@ -250,13 +258,12 @@ function surfaceintegral( end """ - surfaceintegral(f, triangle::Meshes.Triangle, ::GaussKronrod) + integral(f, triangle::Meshes.Triangle, ::GaussKronrod) -Like [`surfaceintegral`](@ref) but integrates over the surface of a `triangle` -using a Gauss-Kronrod quadrature rule along each barycentric dimension of the -triangle. +Like [`integral`](@ref) but integrates over the surface of a `triangle` using nested +Gauss-Kronrod quadrature rules along each barycentric dimension of the triangle. """ -function surfaceintegral( +function integral( f, triangle::Meshes.Ngon{3,Dim,T}, settings::GaussKronrod @@ -273,7 +280,7 @@ function surfaceintegral( return 2.0 * area(triangle) .* outerintegral end -function surfaceintegral( +function integral( f, sphere::Meshes.Sphere{3,T}, settings::GaussKronrod @@ -294,7 +301,7 @@ end # HCubature ################################################################################ -function surfaceintegral( +function integral( f, disk::Meshes.Ball{2,T}, settings::HAdaptiveCubature @@ -314,7 +321,7 @@ function surfaceintegral( return (2π * disk.radius) .* intval end -function surfaceintegral( +function integral( f, box::Meshes.Box{2,T}, settings::HAdaptiveCubature @@ -331,7 +338,7 @@ function surfaceintegral( return area(box) .* intval end -function surfaceintegral( +function integral( f, disk::Meshes.Disk{T}, settings::HAdaptiveCubature @@ -353,13 +360,13 @@ function surfaceintegral( end """ - surfaceintegral(f, triangle::Meshes.Triangle, ::GaussKronrod) + integral(f, triangle::Meshes.Triangle, ::GaussKronrod) -Like [`surfaceintegral`](@ref) but integrates over the surface of a `triangle` -by transforming the triangle into a polar-barycentric coordinate system and -using an h-adaptive cubature rule. +Like [`integral`](@ref) but integrates over the surface of a `triangle` by +transforming the triangle into a polar-barycentric coordinate system and using +an h-adaptive cubature rule. """ -function surfaceintegral( +function integral( f, triangle::Meshes.Ngon{3,Dim,T}, settings::HAdaptiveCubature @@ -385,7 +392,7 @@ function surfaceintegral( return 2.0 * area(triangle) .* intval end -function surfaceintegral( +function integral( f, sphere::Meshes.Sphere{3,T}, settings::HAdaptiveCubature From 39659882afc96782825a8ae70f17ff2ffb401a87 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 25 Feb 2024 11:48:05 -0500 Subject: [PATCH 03/10] Implement renamed functions --- src/integral_line.jl | 67 +++++++++++++++++++++++------------------ src/integral_surface.jl | 8 ----- 2 files changed, 38 insertions(+), 37 deletions(-) diff --git a/src/integral_line.jl b/src/integral_line.jl index b58882a4..47a80c16 100644 --- a/src/integral_line.jl +++ b/src/integral_line.jl @@ -2,7 +2,7 @@ # Common Methods ################################################################################ -function lineintegral( +function integral( f::F, ring::Meshes.Ring, settings::I @@ -11,7 +11,7 @@ function lineintegral( return sum(segment -> lineintegral(f, segment, settings), segments(ring)) end -function lineintegral( +function integral( f::F, rope::Meshes.Rope, settings::I @@ -20,22 +20,31 @@ function lineintegral( return sum(segment -> lineintegral(f, segment, settings), segments(rope)) end +function lineintegral( + f::F, + curve::Meshes.BezierCurve{Dim,T,V}, + settings::I; + alg::Meshes.BezierEvalMethod=Meshes.Horner() +) where {F<:Function, Dim, T, V, I<:IntegrationAlgorithm} + return integral(f, curve, settings; alg=alg) +end + ################################################################################ # Gauss-Legendre ################################################################################ """ - lineintegral(f, curve::Meshes.BezierCurve, ::GaussLegendre; alg=Meshes.Horner()) + integral(f, curve::Meshes.BezierCurve, ::GaussLegendre; alg=Meshes.Horner()) -Like [`lineintegral`](@ref) but integrates along the domain defined a `curve`. -By default this uses Horner's method to improve performance when parameterizing +Like [`integral`](@ref) but integrates along the domain defined a `curve`. By +default this uses Horner's method to improve performance when parameterizing the `curve` at the expense of a small loss of precision. Additional accuracy can be obtained by specifying the use of DeCasteljau's algorithm instead with `alg=Meshes.DeCasteljau()` but can come at a steep cost in memory allocations, especially for curves with a large number of control points. """ -function lineintegral( +function integral( f::F, curve::Meshes.BezierCurve{Dim,T,V}, settings::GaussLegendre; @@ -55,7 +64,7 @@ function lineintegral( return 0.5 * length(curve) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs)) end -function lineintegral( +function integral( f::F, line::Meshes.Box{1,T}, settings::GaussLegendre @@ -75,7 +84,7 @@ function lineintegral( return 0.5 * length(line) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs)) end -function lineintegral( +function integral( f::F, circle::Meshes.Circle{T}, settings::GaussLegendre @@ -96,7 +105,7 @@ function lineintegral( return 0.5 * length(circle) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs)) end -function lineintegral( +function integral( f::F, line::Meshes.Line{Dim,T}, settings::GaussLegendre @@ -118,7 +127,7 @@ function lineintegral( return sum(w .* integrand(x) for (w,x) in zip(ws, xs)) end -function lineintegral( +function integral( f::F, segment::Meshes.Segment{Dim,T}, settings::GaussLegendre @@ -137,7 +146,7 @@ function lineintegral( return 0.5 * length(segment) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs)) end -function lineintegral( +function integral( f::F, circle::Meshes.Sphere{2,T}, settings::GaussLegendre @@ -164,16 +173,16 @@ end ################################################################################ """ - lineintegral(f, curve::BezierCurve, ::GaussKronrod; alg=Horner(), kws...) + integral(f, curve::BezierCurve, ::GaussKronrod; alg=Horner(), kws...) -Like [`lineintegral`](@ref) but integrates along the domain defined a `curve`. -By default this uses Horner's method to improve performance when parameterizing +Like [`integral`](@ref) but integrates along the domain defined a `curve`. By +default this uses Horner's method to improve performance when parameterizing the `curve` at the expense of a small loss of precision. Additional accuracy can be obtained by specifying the use of DeCasteljau's algorithm instead with `alg=Meshes.DeCasteljau()` but can come at a steep cost in memory allocations, especially for curves with a large number of control points. """ -function lineintegral( +function integral( f::F, curve::Meshes.BezierCurve{Dim,T,V}, settings::GaussKronrod; @@ -187,7 +196,7 @@ function lineintegral( return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, line::Meshes.Box{1,T}, settings::GaussKronrod @@ -201,7 +210,7 @@ function lineintegral( return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, circle::Meshes.Circle{T}, settings::GaussKronrod @@ -215,7 +224,7 @@ function lineintegral( return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, line::Meshes.Line{Dim,T}, settings::GaussKronrod @@ -231,7 +240,7 @@ function lineintegral( return QuadGK.quadgk(t -> f(point(t)), -Inf, Inf; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, segment::Meshes.Segment{Dim,T}, settings::GaussKronrod @@ -244,7 +253,7 @@ function lineintegral( return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, circle::Meshes.Sphere{2,T}, settings::GaussKronrod @@ -264,16 +273,16 @@ end ################################################################################ """ - lineintegral(f, curve::BezierCurve, ::HAdaptiveCubature; alg=Horner(), kws...) + integral(f, curve::BezierCurve, ::HAdaptiveCubature; alg=Horner(), kws...) -Like [`lineintegral`](@ref) but integrates along the domain defined a `curve`. -By default this uses Horner's method to improve performance when parameterizing +Like [`integral`](@ref) but integrates along the domain defined a `curve`. By +default this uses Horner's method to improve performance when parameterizing the `curve` at the expense of a small loss of precision. Additional accuracy can be obtained by specifying the use of DeCasteljau's algorithm instead with `alg=Meshes.DeCasteljau()` but can come at a steep cost in memory allocations, especially for curves with a large number of control points. """ -function lineintegral( +function integral( f::F, curve::Meshes.BezierCurve{Dim,T,V}, settings::HAdaptiveCubature; @@ -287,7 +296,7 @@ function lineintegral( return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, line::Meshes.Box{1,T}, settings::HAdaptiveCubature @@ -301,7 +310,7 @@ function lineintegral( return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, circle::Meshes.Circle{T}, settings::HAdaptiveCubature @@ -315,7 +324,7 @@ function lineintegral( return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, line::Meshes.Line{Dim,T}, settings::HAdaptiveCubature @@ -335,7 +344,7 @@ function lineintegral( return hcubature(integrand, [-1], [1]; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, segment::Meshes.Segment{Dim,T}, settings::HAdaptiveCubature @@ -348,7 +357,7 @@ function lineintegral( return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1] end -function lineintegral( +function integral( f::F, circle::Meshes.Sphere{2,T}, settings::HAdaptiveCubature diff --git a/src/integral_surface.jl b/src/integral_surface.jl index fa034146..374e396c 100644 --- a/src/integral_surface.jl +++ b/src/integral_surface.jl @@ -1,11 +1,3 @@ -function lineintegral(f, geometry, settings) - dim = paramdim(geometry) - if == 1 - return integral(f, geometry, settings) - else - error("Unable to perform line integral on a geometry with $dim parametric dimensions.") -end - ################################################################################ # Gauss-Legendre ################################################################################ From 47a9456fa14e4a7d4c3c5917d5a5d03c0de01276 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 25 Feb 2024 11:49:34 -0500 Subject: [PATCH 04/10] Bump version and update test matrix --- Project.toml | 2 +- test/runtests.jl | 32 ++++++++++++++++---------------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index 236c427c..8204973c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MeshIntegrals" uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47" authors = ["Mike Ingold "] -version = "0.9.6" +version = "0.10.0" [deps] FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" diff --git a/test/runtests.jl b/test/runtests.jl index e75bdf5b..73b4a7d8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -98,27 +98,27 @@ end SUPPORT_MATRIX = [ # Name, example, integral,line,surface,volume, GaussLegendre,GaussKronrod,HAdaptiveCubature - SupportItem("BezierCurve", bezier, 0, 1, 0, 0, 1, 1, 1), - SupportItem("Box{1,T}", box1d, 0, 1, 0, 0, 1, 1, 1), - SupportItem("Circle", circle, 0, 1, 0, 0, 1, 1, 1), + SupportItem("BezierCurve", bezier, 1, 1, 0, 0, 1, 1, 1), + SupportItem("Box{1,T}", box1d, 1, 1, 0, 0, 1, 1, 1), + SupportItem("Circle", circle, 1, 1, 0, 0, 1, 1, 1), # Line -- custom test - SupportItem("Ring", ring_rect, 0, 1, 0, 0, 1, 1, 1), - SupportItem("Rope", rope_rect, 0, 1, 0, 0, 1, 1, 1), - SupportItem("Segment", seg_ne, 0, 1, 0, 0, 1, 1, 1), - SupportItem("Sphere{2,T}", sphere2d, 0, 1, 0, 0, 1, 1, 1), - - SupportItem("Ball{2,T}", ball2d, 0, 0, 1, 0, 1, 1, 1), - SupportItem("Box{2,T}", box2d, 0, 0, 1, 0, 1, 1, 1), - SupportItem("CylinderSurface", cylsurf, 0, 0, 1, 0, 0, 1, 0), - SupportItem("Disk", disk, 0, 0, 1, 0, 1, 1, 1), + SupportItem("Ring", ring_rect, 1, 1, 0, 0, 1, 1, 1), + SupportItem("Rope", rope_rect, 1, 1, 0, 0, 1, 1, 1), + SupportItem("Segment", seg_ne, 1, 1, 0, 0, 1, 1, 1), + SupportItem("Sphere{2,T}", sphere2d, 1, 1, 0, 0, 1, 1, 1), + + SupportItem("Ball{2,T}", ball2d, 1, 0, 1, 0, 1, 1, 1), + SupportItem("Box{2,T}", box2d, 1, 0, 1, 0, 1, 1, 1), + SupportItem("CylinderSurface", cylsurf, 1, 0, 1, 0, 0, 1, 0), + SupportItem("Disk", disk, 1, 0, 1, 0, 1, 1, 1), # ParaboloidSurface -- not yet supported - SupportItem("Sphere{3,T}", sphere3d, 0, 0, 1, 0, 1, 1, 1), - SupportItem("Triangle", triangle, 0, 0, 1, 0, 1, 1, 1), + SupportItem("Sphere{3,T}", sphere3d, 1, 0, 1, 0, 1, 1, 1), + SupportItem("Triangle", triangle, 1, 0, 1, 0, 1, 1, 1), # Torus -- not yet supported # SimpleMesh -- not yet supported - SupportItem("Ball{3,T}", ball3d, 0, 0, 0, 1, 1, 0, 1), - SupportItem("Box{3,T}", box3d, 0, 0, 0, 1, 1, 0, 1) + SupportItem("Ball{3,T}", ball3d, 1, 0, 0, 1, 1, 0, 1), + SupportItem("Box{3,T}", box3d, 1, 0, 0, 1, 1, 0, 1) ] # Run all integral tests From 8ab9edb505c4faf3759b5b975acbba20786b661b Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 25 Feb 2024 21:16:43 -0500 Subject: [PATCH 05/10] Update naming --- src/integral_volume.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/integral_volume.jl b/src/integral_volume.jl index e62ea9be..b337fd75 100644 --- a/src/integral_volume.jl +++ b/src/integral_volume.jl @@ -2,7 +2,7 @@ # Gauss-Legendre ################################################################################ -function volumeintegral( +function integral( f, ball::Meshes.Ball{3,T}, settings::GaussLegendre @@ -27,7 +27,7 @@ function volumeintegral( return (1/8) * 2π^2 * R^3 .* sum(g, zip(wws,xxs)) end -function volumeintegral( +function integral( f, box::Meshes.Box{3,T}, settings::GaussLegendre @@ -59,7 +59,7 @@ end # HCubature ################################################################################ -function volumeintegral( +function integral( f, ball::Meshes.Ball{3,T}, settings::HAdaptiveCubature @@ -76,7 +76,7 @@ function volumeintegral( return 2π^2 * R^3 .* intval end -function volumeintegral( +function integral( f, box::Meshes.Box{3,T}, settings::HAdaptiveCubature From 4a37150d9cc65577b7025dfc74c4e7d5dd90cfa5 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 25 Feb 2024 21:27:48 -0500 Subject: [PATCH 06/10] Bugfix --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 73b4a7d8..26796e0c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,7 +31,7 @@ function integraltest(intf, geometry, rule, supported) @test intf(f, geometry, rule) ≈ measure(geometry) @test intf(fv, geometry, rule) ≈ fill(measure(geometry),3) else - @test_throws MethodError intf(f, geometry, rule) + @test_throws "Unable to perform" intf(f, geometry, rule) end end From fcb9e746e1adaf59d0df297e96dfd3a7e282320e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Mar 2024 08:26:20 -0500 Subject: [PATCH 07/10] Add debugging for test generation --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 26796e0c..7040649e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -55,6 +55,7 @@ function autotest(item::SupportItem) # For each enabled solver type, run the test suite @testset "$(item.name)" begin for ((method,methodsupport), (alg,algsupport)) in itemsupport + @info method, methodsupport, alg, algsupport integraltest(method, item.geometry, alg, methodsupport && algsupport) end end From 5c0ef4383d5d3b76a31f5febaad39c270e968adf Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Mar 2024 08:32:59 -0500 Subject: [PATCH 08/10] More debugging --- test/runtests.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7040649e..42486d0a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -55,8 +55,9 @@ function autotest(item::SupportItem) # For each enabled solver type, run the test suite @testset "$(item.name)" begin for ((method,methodsupport), (alg,algsupport)) in itemsupport - @info method, methodsupport, alg, algsupport - integraltest(method, item.geometry, alg, methodsupport && algsupport) + support = methodsupport && algsupport + @info item.name, method, methodsupport, alg, algsupport, support + integraltest(method, item.geometry, alg, support) end end end From e8959412f75f07e833ef6217c71a7d411d40c80f Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Mar 2024 09:05:39 -0500 Subject: [PATCH 09/10] Remove debugging and update messaging for not-supported methods --- src/integral.jl | 6 +++--- src/integral_surface.jl | 18 ++++++++++++++++++ src/integral_volume.jl | 21 +++++++++++++++++++++ test/runtests.jl | 6 ++---- 4 files changed, 44 insertions(+), 7 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index ea0a18dc..24af9fbf 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -84,7 +84,7 @@ function lineintegral( if dim == 1 return integral(f, geometry, settings) else - error("Unable to perform line integral on a geometry with $dim parametric dimensions.") + error("Performing a line integral on a geometry with $dim parametric dimensions not supported.") end end @@ -103,7 +103,7 @@ function surfaceintegral( if dim == 2 return integral(f, geometry, settings) else - error("Unable to perform surface integral on a geometry with $dim parametric dimensions.") + error("Performing a surface integral on a geometry with $dim parametric dimensions not supported.") end end @@ -123,6 +123,6 @@ function volumeintegral( if dim == 3 return integral(f, geometry, settings) else - error("Unable to perform surface integral on a geometry with $dim parametric dimensions.") + error("Performing a volume integral on a geometry with $dim parametric dimensions not supported.") end end diff --git a/src/integral_surface.jl b/src/integral_surface.jl index 374e396c..00365383 100644 --- a/src/integral_surface.jl +++ b/src/integral_surface.jl @@ -54,6 +54,15 @@ function integral( return 0.25 * area(box) .* sum(g, zip(wws,xxs)) end +function integral( + f, + cyl::Meshes.CylinderSurface{T}, + settings::GaussLegendre +) where {T} + error("Integrating a CylinderSurface{T} with GaussLegendre not supported.") + # Planned to support in the future +end + function integral( f, disk::Meshes.Disk{T}, @@ -330,6 +339,15 @@ function integral( return area(box) .* intval end +function integral( + f, + cyl::Meshes.CylinderSurface{T}, + settings::HAdaptiveCubature +) where {T} + error("Integrating a CylinderSurface{T} with HAdaptiveCubature not supported.") + # Planned to support in the future +end + function integral( f, disk::Meshes.Disk{T}, diff --git a/src/integral_volume.jl b/src/integral_volume.jl index b337fd75..ac219435 100644 --- a/src/integral_volume.jl +++ b/src/integral_volume.jl @@ -55,6 +55,27 @@ function integral( end +################################################################################ +# GaussKronrod +################################################################################ + +function integral( + f, + ball::Meshes.Ball{3,T}, + settings::GaussKronrod +) where {T} + error("Integrating a Ball{3,T} with GaussKronrod not supported.") +end + +function integral( + f, + box::Meshes.Box{3,T}, + settings::GaussKronrod +) where {T} + error("Integrating a Ball{3,T} with GaussKronrod not supported.") +end + + ################################################################################ # HCubature ################################################################################ diff --git a/test/runtests.jl b/test/runtests.jl index 42486d0a..0277d01c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,7 +31,7 @@ function integraltest(intf, geometry, rule, supported) @test intf(f, geometry, rule) ≈ measure(geometry) @test intf(fv, geometry, rule) ≈ fill(measure(geometry),3) else - @test_throws "Unable to perform" intf(f, geometry, rule) + @test_throws "not supported" intf(f, geometry, rule) end end @@ -55,9 +55,7 @@ function autotest(item::SupportItem) # For each enabled solver type, run the test suite @testset "$(item.name)" begin for ((method,methodsupport), (alg,algsupport)) in itemsupport - support = methodsupport && algsupport - @info item.name, method, methodsupport, alg, algsupport, support - integraltest(method, item.geometry, alg, support) + integraltest(method, item.geometry, alg, methodsupport && algsupport) end end end From 3943a0056e3c9d0fee7adc24b78fbab12bd5e53d Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Mar 2024 09:17:35 -0500 Subject: [PATCH 10/10] Update README --- README.md | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index ec97344c..14205f74 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,14 @@ # MeshIntegrals.jl -This package implements methods for computing integrals over geometric polytopes +This package implements methods for numerically-computing integrals over geometric polytopes from [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) using: -- Gauss-Legendre quadrature rules from [**FastGaussQuadrature.jl**](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) -- H-adaptive Gauss-Kronrod quadrature rules from [**QuadGK.jl**](https://github.com/JuliaMath/QuadGK.jl) -- H-adaptive cubature rules from [**HCubature.jl**](https://github.com/JuliaMath/HCubature.jl) +- Gauss-Legendre quadrature rules from [**FastGaussQuadrature.jl**](https://github.com/JuliaApproximation/FastGaussQuadrature.jl): `GaussLegendre(n)` +- H-adaptive Gauss-Kronrod quadrature rules from [**QuadGK.jl**](https://github.com/JuliaMath/QuadGK.jl): `GaussKronrod(kwargs...)` +- H-adaptive cubature rules from [**HCubature.jl**](https://github.com/JuliaMath/HCubature.jl): `HAdaptiveCubature(kwargs...)` + +Functions available: +- `integral(f, geometry, ::IntegrationAlgorithm)`: integrates a function `f` over a domain defined by `geometry` using a particular +- `lineintegral`, `surfaceintegral`, and `volumeintegral` are available as aliases for `integral` that first verify that `geometry` has the appropriate number of parametric dimensions Methods are tested to ensure compatibility with - Meshes.jl geometries with **Unitful.jl** coordinate types, e.g. `Point(1.0u"m", 2.0u"m")`