diff --git a/Project.toml b/Project.toml index 0cfb8258..5addc43a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,6 @@ name = "MeshIntegrals" uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47" -authors = ["Mike Ingold "] -version = "0.13.4" +version = "0.13.5" [deps] CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb" diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 74dd4925..48825a71 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -11,11 +11,23 @@ module MeshIntegrals include("utils.jl") export jacobian, derivative, unitdirection + include("integration_algorithms.jl") + export GaussKronrod, GaussLegendre, HAdaptiveCubature + include("integral.jl") include("integral_aliases.jl") - include("integral_line.jl") - include("integral_surface.jl") - include("integral_volume.jl") - export GaussKronrod, GaussLegendre, HAdaptiveCubature export integral, lineintegral, surfaceintegral, volumeintegral + + # Integration methods specialized for particular geometries + include("specializations/BezierCurve.jl") + include("specializations/ConeSurface.jl") + include("specializations/CylinderSurface.jl") + include("specializations/FrustumSurface.jl") + include("specializations/Line.jl") + include("specializations/Plane.jl") + include("specializations/Ray.jl") + include("specializations/Ring.jl") + include("specializations/Rope.jl") + include("specializations/Tetrahedron.jl") + include("specializations/Triangle.jl") end diff --git a/src/integral.jl b/src/integral.jl index c475fd80..f170164f 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -1,48 +1,3 @@ -################################################################################ -# Integration Algorithms -################################################################################ - -abstract type IntegrationAlgorithm end - -""" - GaussKronrod(kwargs...) - -Numerically integrate using the h-adaptive Gauss-Kronrod quadrature rule implemented -by QuadGK.jl. All standard `QuadGK.quadgk` keyword arguments are supported. -""" -struct GaussKronrod <: IntegrationAlgorithm - kwargs - GaussKronrod(; kwargs...) = new(kwargs) -end - -""" - GaussLegendre(n) - -Numerically integrate using an `n`'th-order Gauss-Legendre quadrature rule. nodes -and weights are efficiently calculated using FastGaussQuadrature.jl. - -So long as the integrand function can be well-approximated by a polynomial of -order `2n-1`, this method should yield results with 16-digit accuracy in `O(n)` -time. If the function is know to have some periodic content, then `n` should -(at a minimum) be greater than the expected number of periods over the geometry, -e.g. `length(geometry)/lambda`. -""" -struct GaussLegendre <: IntegrationAlgorithm - n::Int64 -end - -""" - GaussKronrod(kwargs...) - -Numerically integrate areas and surfaces using the h-adaptive cubature rule -implemented by HCubature.jl. All standard `HCubature.hcubature` keyword -arguments are supported. -""" -struct HAdaptiveCubature <: IntegrationAlgorithm - kwargs - HAdaptiveCubature(; kwargs...) = new(kwargs) -end - ################################################################################ # Master Integral Function ################################################################################ @@ -68,57 +23,38 @@ contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision """ function integral end -# If only f and geometry are specified, default to HAdaptiveCubature +# If only f and geometry are specified, select default algorithm function integral( f::F, geometry::G ) where {F<:Function, G<:Meshes.Geometry} - _integral(f, geometry, HAdaptiveCubature()) -end - -# If algorithm is HAdaptiveCubature, use the generalized n-dim solution -function integral( - f::F, - geometry::G, - settings::HAdaptiveCubature -) where {F<:Function, G<:Meshes.Geometry} - _integral(f, geometry, settings) + N = Meshes.paramdim(geometry) + rule = (N == 1) ? GaussKronrod() : HAdaptiveCubature() + _integral(f, geometry, rule) end -# If algorithm is HAdaptiveCubature, and FP specified, use the generalized n-dim solution +# with algorithm and T specified function integral( f::F, geometry::G, - settings::HAdaptiveCubature, - FP::Type{T} -) where {F<:Function, G<:Meshes.Geometry, T<:AbstractFloat} + settings::I, + FP::Type{T} = Float64 +) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat} _integral(f, geometry, settings, FP) end -# If algorithm is not HAdaptiveCubature, specialize on number of dimensions -function integral( - f::F, - geometry::G, - settings::I -) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm} - # Run the appropriate integral type - Dim = Meshes.paramdim(geometry) - if Dim == 1 - return _integral_1d(f, geometry, settings) - elseif Dim == 2 - return _integral_2d(f, geometry, settings) - elseif Dim == 3 - return _integral_3d(f, geometry, settings) - end -end -# If algorithm is not HAdaptiveCubature, and FP specified, specialize on number of dimensions -function integral( - f::F, - geometry::G, - settings::I, - FP::Type{T} -) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat} +################################################################################ +# Generalized (n-Dimensional) Worker Methods +################################################################################ + +# GaussKronrod +function _integral( + f, + geometry, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {T<:AbstractFloat} # Run the appropriate integral type Dim = Meshes.paramdim(geometry) if Dim == 1 @@ -130,12 +66,32 @@ function integral( end end +# GaussLegendre +function _integral( + f, + geometry, + settings::GaussLegendre, + FP::Type{T} = Float64 +) where {T<:AbstractFloat} + N = Meshes.paramdim(geometry) + + # Get Gauss-Legendre nodes and weights for a region [-1,1]^N + xs, ws = _gausslegendre(FP, settings.n) + weights = Iterators.product(ntuple(Returns(ws), N)...) + nodes = Iterators.product(ntuple(Returns(xs), N)...) -################################################################################ -# Generalized (n-Dimensional) Worker Methods -################################################################################ + # Domain transformation: x [-1,1] ↦ t [0,1] + t(x) = FP(1//2) * x + FP(1//2) -# General solution for HAdaptiveCubature methods + function integrand((weights, nodes)) + ts = t.(nodes) + prod(weights) * f(geometry(ts...)) * differential(geometry, ts) + end + + return FP(1//(2^N)) .* sum(integrand, zip(weights, nodes)) +end + +# HAdaptiveCubature function _integral( f, geometry, @@ -158,3 +114,38 @@ function _integral( # Reapply units return value .* integrandunits end + +################################################################################ +# Specialized GaussKronrod Methods +################################################################################ + +function _integral_1d( + f, + geometry, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {T<:AbstractFloat} + integrand(t) = f(geometry(t)) * differential(geometry, [t]) + return QuadGK.quadgk(integrand, FP(0), FP(1); settings.kwargs...)[1] +end + +function _integral_2d( + f, + geometry2d, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {T<:AbstractFloat} + integrand(u,v) = f(geometry2d(u,v)) * differential(geometry2d, [u,v]) + ∫₁(v) = QuadGK.quadgk(u -> integrand(u,v), FP(0), FP(1); settings.kwargs...)[1] + return QuadGK.quadgk(v -> ∫₁(v), FP(0), FP(1); settings.kwargs...)[1] +end + +# Integrating volumes with GaussKronrod not supported by default +function _integral_3d( + f, + geometry, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {T<:AbstractFloat} + error("Integrating this volume type with GaussKronrod not supported.") +end diff --git a/src/integral_line.jl b/src/integral_line.jl deleted file mode 100644 index 39939b9a..00000000 --- a/src/integral_line.jl +++ /dev/null @@ -1,381 +0,0 @@ -################################################################################ -# Generalized 1D Methods -################################################################################ - -function _integral_1d( - f, - geometry, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {T<:AbstractFloat} - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, settings.n) - - # Change of variables: x [-1,1] ↦ t [0,1] - t(x) = FP(1/2) * x + FP(1/2) - - # Integrate f along the geometry and apply a domain-correction factor for [-1,1] ↦ [0, 1] - integrand((w,x)) = w * f(geometry(t(x))) * differential(geometry, [t(x)]) - return FP(1/2) * sum(integrand, zip(ws, xs)) -end - -function _integral_1d( - f, - geometry, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {T<:AbstractFloat} - integrand(t) = f(geometry(t)) * differential(geometry, [t]) - return QuadGK.quadgk(integrand, FP(0), FP(1); settings.kwargs...)[1] -end - -function _integral_1d( - f, - geometry, - settings::HAdaptiveCubature, -) - return _integral(f, geometry, settings) -end - -function _integral_1d( - f, - geometry, - settings::HAdaptiveCubature, - FP::Type{T} -) where {T<:AbstractFloat} - return _integral(f, geometry, settings, FP) -end - - -################################################################################ -# Specialized Methods for BezierCurve -################################################################################ - -function lineintegral( - f::F, - curve::Meshes.BezierCurve, - settings::I, - FP::Type{T}; - alg::Meshes.BezierEvalMethod=Meshes.Horner() -) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} - return integral(f, curve, settings, FP; alg=alg) -end - -function lineintegral( - f::F, - curve::Meshes.BezierCurve, - settings::I; - alg::Meshes.BezierEvalMethod=Meshes.Horner() -) where {F<:Function, I<:IntegrationAlgorithm} - return integral(f, curve, settings; alg=alg) -end - -""" - integral(f, curve::Meshes.BezierCurve, ::GaussLegendre; alg=Meshes.Horner()) - -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 integral( - f::F, - curve::Meshes.BezierCurve, - settings::GaussLegendre, - FP::Type{T} = Float64; - alg::Meshes.BezierEvalMethod=Meshes.Horner() -) where {F<:Function, T<:AbstractFloat} - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, settings.n) - - # Change of variables: x [-1,1] ↦ t [0,1] - t(x) = FP(1/2) * x + FP(1/2) - point(x) = curve(t(x), alg) - - # Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length] - return FP(1/2) * length(curve) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs)) -end - -""" - integral(f, curve::BezierCurve, ::GaussKronrod; alg=Horner(), kws...) - -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 integral( - f::F, - curve::Meshes.BezierCurve, - settings::GaussKronrod, - FP::Type{T} = Float64; - alg::Meshes.BezierEvalMethod=Meshes.Horner() -) where {F<:Function, T<:AbstractFloat} - len = length(curve) - point(t) = curve(t, alg) - return QuadGK.quadgk(t -> len * f(point(t)), FP(0), FP(1); settings.kwargs...)[1] -end - -""" - integral(f, curve::BezierCurve, ::HAdaptiveCubature; alg=Horner(), kws...) - -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 integral( - f::F, - curve::Meshes.BezierCurve, - settings::HAdaptiveCubature, - FP::Type{T} = Float64; - alg::Meshes.BezierEvalMethod=Meshes.Horner() -) where {F<:Function, T<:AbstractFloat} - len = length(curve) - point(t) = curve(t, alg) - integrand(t) = len * f(point(t[1])) - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = fill(FP(0.5),3) - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, FP[0], FP[1]; settings.kwargs...)[1] - - # Reapply units - return value .* integrandunits -end - -################################################################################ -# Specialized Methods for Line -################################################################################ - -function integral( - f::F, - line::Meshes.Line, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, settings.n) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Line - domainunits = _units(line(0)) - integrand(x) = f(line(t(x))) * t′(x) * domainunits - return sum(w .* integrand(x) for (w,x) in zip(ws, xs)) -end - -function integral( - f::F, - line::Meshes.Line, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Normalize the Line s.t. line(t) is distance t from origin point - line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Integrate f along the Line - domainunits = _units(line(0)) - integrand(t) = f(line(t)) * domainunits - return QuadGK.quadgk(integrand, FP(-Inf), FP(Inf); settings.kwargs...)[1] -end - -function integral( - f::F, - line::Meshes.Line, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Normalize the Line s.t. line(t) is distance t from origin point - line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Line - domainunits = _units(line(0)) - integrand(x::AbstractVector) = f(line(t(x[1]))) * t′(x[1]) * domainunits - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = FP[0.5] - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, FP[-1], FP[1]; settings.kwargs...)[1] - - # Reapply units - return value .* integrandunits -end - -################################################################################ -# Specialized Methods for Ray -################################################################################ - -function integral( - f::F, - ray::Meshes.Ray, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, settings.n) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Ray(ray.p, Meshes.unormalize(ray.v)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞) - t₁(x) = FP(1/2) * x + FP(1/2) - t₁′(x) = FP(1/2) - t₂(x) = x / (1 - x^2) - t₂′(x) = (1 + x^2) / (1 - x^2)^2 - t = t₂ ∘ t₁ - t′(x) = t₂′(t₁(x)) * t₁′(x) - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(x) = f(ray(t(x))) * t′(x) * domainunits - return sum(w .* integrand(x) for (w,x) in zip(ws, xs)) -end - -function integral( - f::F, - ray::Meshes.Ray, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Ray(ray.p, Meshes.unormalize(ray.v)) - - # Integrate f along the Ray - domainunits = _units(ray(0)) - return QuadGK.quadgk(t -> f(ray(t)) * domainunits, FP(0), FP(Inf); settings.kwargs...)[1] -end - -function integral( - f::F, - ray::Meshes.Ray, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Ray(ray.p, Meshes.unormalize(ray.v)) - - # Domain transformation: x ∈ [0,1] ↦ t ∈ [0,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(x::AbstractVector) = f(ray(t(x[1]))) * t′(x[1]) * domainunits - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = FP[0.5] - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, FP[0], FP[1]; settings.kwargs...)[1] - - # Reapply units - return value .* integrandunits -end - -################################################################################ -# Specialized Methods for Ring, Rope -################################################################################ - -function integral( - f::F, - ring::Meshes.Ring, - settings::HAdaptiveCubature -) where {F<:Function} - # Convert the Ring into Segments, sum the integrals of those - return sum(segment -> _integral(f, segment, settings), segments(ring)) -end - -function integral( - f::F, - ring::Meshes.Ring, - settings::HAdaptiveCubature, - FP::Type{T} -) where {F<:Function, T<:AbstractFloat} - # Convert the Ring into Segments, sum the integrals of those - return sum(segment -> _integral(f, segment, settings, FP), segments(ring)) -end - -function integral( - f::F, - ring::Meshes.Ring, - settings::I -) where {F<:Function, I<:IntegrationAlgorithm} - # Convert the Ring into Segments, sum the integrals of those - return sum(segment -> integral(f, segment, settings), segments(ring)) -end - -function integral( - f::F, - ring::Meshes.Ring, - settings::I, - FP::Type{T} -) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} - # Convert the Ring into Segments, sum the integrals of those - return sum(segment -> integral(f, segment, settings, FP), segments(ring)) -end - - -function integral( - f::F, - rope::Meshes.Rope, - settings::HAdaptiveCubature -) where {F<:Function} - # Convert the Rope into Segments, sum the integrals of those - return sum(segment -> _integral(f, segment, settings), segments(rope)) -end - -function integral( - f::F, - rope::Meshes.Rope, - settings::HAdaptiveCubature, - FP::Type{T} -) where {F<:Function, T<:AbstractFloat} - # Convert the Rope into Segments, sum the integrals of those - return sum(segment -> _integral(f, segment, settings, FP), segments(rope)) -end - -function integral( - f::F, - rope::Meshes.Rope, - settings::I -) where {F<:Function, I<:IntegrationAlgorithm} - # Convert the Rope into Segments, sum the integrals of those - return sum(segment -> integral(f, segment, settings), segments(rope)) -end - -function integral( - f::F, - rope::Meshes.Rope, - settings::I, - FP::Type{T} -) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} - # Convert the Rope into Segments, sum the integrals of those - return sum(segment -> integral(f, segment, settings, FP), segments(rope)) -end diff --git a/src/integral_surface.jl b/src/integral_surface.jl deleted file mode 100644 index 13ade0cb..00000000 --- a/src/integral_surface.jl +++ /dev/null @@ -1,314 +0,0 @@ -################################################################################ -# Generalized 2D Methods -################################################################################ - -function _integral_2d( - f, - geometry2d, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {T<:AbstractFloat} - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² - xs, ws = _gausslegendre(FP, settings.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Domain transformation: x [-1,1] ↦ t [0,1] - t(x) = FP(1/2) * x + FP(1/2) - point(xi, xj) = geometry2d(t(xi), t(xj)) - - # Integrate f over the geometry - integrand(((wi,wj), (xi,xj))) = wi * wj * f(point(xi,xj)) * differential(geometry2d, t.([xi, xj])) - return FP(1/4) .* sum(integrand, zip(wws,xxs)) -end - -function _integral_2d( - f, - geometry2d, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {T<:AbstractFloat} - integrand(u,v) = f(geometry2d(u,v)) * differential(geometry2d, [u,v]) - ∫₁(v) = QuadGK.quadgk(u -> integrand(u,v), FP(0), FP(1); settings.kwargs...)[1] - return QuadGK.quadgk(v -> ∫₁(v), FP(0), FP(1); settings.kwargs...)[1] -end - - -################################################################################ -# Specialized Methods for ConeSurface, CylinderSurface, and FrustumSurface -################################################################################ - -function integral( - f::F, - cyl::Meshes.CylinderSurface, - settings::I, - FP::Type{T} = Float64 -) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} - # The generic method only parameterizes the sides - sides = _integral_2d(f, cyl, settings, FP) - - # Integrate the Disk at the top - disk_top = Meshes.Disk(cyl.top, cyl.radius) - top = _integral_2d(f, disk_top, settings, FP) - - # Integrate the Disk at the bottom - disk_bottom = Meshes.Disk(cyl.bot, cyl.radius) - bottom = _integral_2d(f, disk_bottom, settings, FP) - - return sides + top + bottom -end - -function integral( - f::F, - cyl::Meshes.CylinderSurface, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # The generic method only parameterizes the sides - sides = _integral(f, cyl, settings, FP) - - # Integrate the Disk at the top - disk_top = Meshes.Disk(cyl.top, cyl.radius) - top = _integral(f, disk_top, settings, FP) - - # Integrate the Disk at the bottom - disk_bottom = Meshes.Disk(cyl.bot, cyl.radius) - bottom = _integral(f, disk_bottom, settings, FP) - - return sides + top + bottom -end - -function integral( - f::F, - cone::Meshes.ConeSurface, - settings::I, - FP::Type{T} = Float64 -) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} - # The generic method only parameterizes the sides - sides = _integral_2d(f, cone, settings, FP) - - # Integrate the Disk at the base - base = _integral_2d(f, cone.base, settings, FP) - - return sides + base -end - -function integral( - f::F, - cone::Meshes.ConeSurface, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # The generic method only parameterizes the sides - sides = _integral(f, cone, settings, FP) - - # Integrate the Disk at the base - base = _integral(f, cone.base, settings, FP) - - return sides + base -end - -function integral( - f::F, - frust::Meshes.FrustumSurface, - settings::I, - FP::Type{T} = Float64 -) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} - # The generic method only parameterizes the sides - sides = _integral_2d(f, frust, settings, FP) - - # Integrate the Disks at the top and bottom - top = _integral_2d(f, frust.top, settings, FP) - bottom = _integral_2d(f, frust.bot, settings, FP) - - return sides + top + bottom -end - -function integral( - f::F, - frust::Meshes.FrustumSurface, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # The generic method only parameterizes the sides - sides = _integral(f, frust, settings, FP) - - # Integrate the Disks at the top and bottom - top = _integral(f, frust.top, settings, FP) - bottom = _integral(f, frust.bot, settings, FP) - - return sides + top + bottom -end - -################################################################################ -# Specialized Methods for Plane -################################################################################ - -function integral( - f::F, - plane::Meshes.Plane, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² - xs, ws = _gausslegendre(FP, settings.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Normalize the Plane's orthogonal vectors - plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f over the Plane - domainunits = _units(plane(0,0)) - integrand(((wi,wj), (xi,xj))) = wi * wj * f(plane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2 - return sum(integrand, zip(wws,xxs)) -end - -function integral( - f::F, - plane::Meshes.Plane, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Normalize the Plane's orthogonal vectors - plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v)) - - # Integrate f over the Plane - domainunits = _units(plane(0,0)) - inner∫(v) = QuadGK.quadgk(u -> f(plane(u,v)) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1] - return QuadGK.quadgk(v -> inner∫(v) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1] -end - -function integral( - f::F, - plane::Meshes.Plane, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Normalize the Plane's orthogonal vectors - plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f over the Plane - domainunits = _units(plane(0,0)) - integrand(x::AbstractVector) = f(plane(t(x[1]), t(x[2]))) * t′(x[1]) * t′(x[2]) * domainunits^2 - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = FP[0.5, 0.5] - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, FP[-1,-1], FP[1,1]; settings.kwargs...)[1] - - # Reapply units - return value .* integrandunits -end - -################################################################################ -# Specialized Methods for Triangle -################################################################################ - -""" - integral(f, triangle::Meshes.Triangle, ::GaussLegendre) - -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 integral( - f::F, - triangle::Meshes.Ngon{3}, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 - xs, ws = _gausslegendre(FP, settings.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Domain transformations: - # xᵢ [-1,1] ↦ R [0,1] - # xⱼ [-1,1] ↦ φ [0,π/2] - uR(xᵢ) = T(1/2) * (xᵢ + 1) - uφ(xⱼ) = T(π/4) * (xⱼ + 1) - - # Integrate the Barycentric triangle by transforming it into polar coordinates - # with a modified radius - # R = r ( sinφ + cosφ ) - # s.t. integration bounds become rectangular - # R ∈ [0, 1] and φ ∈ [0, π/2] - function integrand(((wᵢ,wⱼ), (xᵢ,xⱼ))) - R = uR(xᵢ) - φ = uφ(xⱼ) - a,b = sincos(φ) - u = R * (1 - a / (a + b)) - v = R * (1 - b / (a + b)) - return wᵢ * wⱼ * f(triangle(u, v)) * R / (a + b)^2 - end - - # Calculate 2D Gauss-Legendre integral over modified-polar-Barycentric coordinates - # Apply a linear domain-correction factor - return FP(π/4) * area(triangle) .* sum(integrand, zip(wws,xxs)) -end - -""" - integral(f, triangle::Meshes.Triangle, ::GaussKronrod) - -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 integral( - f::F, - triangle::Meshes.Ngon{3}, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0) - # i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du - inner∫(u) = QuadGK.quadgk(v -> f(triangle(u,v)), FP(0), FP(1-u); settings.kwargs...)[1] - outer∫ = QuadGK.quadgk(inner∫, FP(0), FP(1); settings.kwargs...)[1] - - # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * area(triangle) .* outer∫ -end - -""" - integral(f, triangle::Meshes.Triangle, ::GaussKronrod) - -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 integral( - f::F, - triangle::Meshes.Ngon{3}, - settings::HAdaptiveCubature, - FP::Type{T} = Float64 -) where {F<:Function, T<:AbstractFloat} - # Integrate the Barycentric triangle by transforming it into polar coordinates - # with a modified radius - # R = r ( sinφ + cosφ ) - # s.t. integration bounds become rectangular - # R ∈ [0, 1] and φ ∈ [0, π/2] - function integrand(Rφ) - R,φ = Rφ - a,b = sincos(φ) - u = R * (1 - a / (a + b)) - v = R * (1 - b / (a + b)) - return f(triangle(u, v)) * R / (a + b)^2 - end - intval = HCubature.hcubature(integrand, FP[0, 0], FP[1, π/2], settings.kwargs...)[1] - - # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * area(triangle) .* intval -end diff --git a/src/integration_algorithms.jl b/src/integration_algorithms.jl new file mode 100644 index 00000000..8a8eb0b5 --- /dev/null +++ b/src/integration_algorithms.jl @@ -0,0 +1,44 @@ +################################################################################ +# Integration Algorithms +################################################################################ + +abstract type IntegrationAlgorithm end + +""" + GaussKronrod(kwargs...) + +Numerically integrate using the h-adaptive Gauss-Kronrod quadrature rule implemented +by QuadGK.jl. All standard `QuadGK.quadgk` keyword arguments are supported. +""" +struct GaussKronrod <: IntegrationAlgorithm + kwargs + GaussKronrod(; kwargs...) = new(kwargs) +end + +""" + GaussLegendre(n) + +Numerically integrate using an `n`'th-order Gauss-Legendre quadrature rule. Nodes +and weights are efficiently calculated using FastGaussQuadrature.jl. + +So long as the integrand function can be well-approximated by a polynomial of +order `2n-1`, this method should yield results with 16-digit accuracy in `O(n)` +time. If the function is know to have some periodic content, then `n` should +(at a minimum) be greater than the expected number of periods over the geometry, +e.g. `length(geometry)/lambda`. +""" +struct GaussLegendre <: IntegrationAlgorithm + n::Int64 +end + +""" + HAdaptiveCubature(kwargs...) + +Numerically integrate areas and surfaces using the h-adaptive cubature rule +implemented by HCubature.jl. All standard `HCubature.hcubature` keyword +arguments are supported. +""" +struct HAdaptiveCubature <: IntegrationAlgorithm + kwargs + HAdaptiveCubature(; kwargs...) = new(kwargs) +end diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl new file mode 100644 index 00000000..f85b9116 --- /dev/null +++ b/src/specializations/BezierCurve.jl @@ -0,0 +1,106 @@ +################################################################################ +# Specialized Methods for BezierCurve +################################################################################ + +function lineintegral( + f::F, + curve::Meshes.BezierCurve, + settings::I, + FP::Type{T}; + alg::Meshes.BezierEvalMethod=Meshes.Horner() +) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} + return integral(f, curve, settings, FP; alg=alg) +end + +function lineintegral( + f::F, + curve::Meshes.BezierCurve, + settings::I; + alg::Meshes.BezierEvalMethod=Meshes.Horner() +) where {F<:Function, I<:IntegrationAlgorithm} + return integral(f, curve, settings; alg=alg) +end + +""" + integral(f, curve::Meshes.BezierCurve, ::GaussLegendre; alg=Meshes.Horner()) + +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 integral( + f::F, + curve::Meshes.BezierCurve, + settings::GaussLegendre, + FP::Type{T} = Float64; + alg::Meshes.BezierEvalMethod=Meshes.Horner() +) where {F<:Function, T<:AbstractFloat} + # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] + xs, ws = _gausslegendre(FP, settings.n) + + # Change of variables: x [-1,1] ↦ t [0,1] + t(x) = FP(1/2) * x + FP(1/2) + point(x) = curve(t(x), alg) + + # Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length] + return FP(1/2) * length(curve) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs)) +end + +""" + integral(f, curve::BezierCurve, ::GaussKronrod; alg=Horner(), kws...) + +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 integral( + f::F, + curve::Meshes.BezierCurve, + settings::GaussKronrod, + FP::Type{T} = Float64; + alg::Meshes.BezierEvalMethod=Meshes.Horner() +) where {F<:Function, T<:AbstractFloat} + len = length(curve) + point(t) = curve(t, alg) + return QuadGK.quadgk(t -> len * f(point(t)), FP(0), FP(1); settings.kwargs...)[1] +end + +""" + integral(f, curve::BezierCurve, ::HAdaptiveCubature; alg=Horner(), kws...) + +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 integral( + f::F, + curve::Meshes.BezierCurve, + settings::HAdaptiveCubature, + FP::Type{T} = Float64; + alg::Meshes.BezierEvalMethod=Meshes.Horner() +) where {F<:Function, T<:AbstractFloat} + len = length(curve) + point(t) = curve(t, alg) + integrand(t) = len * f(point(t[1])) + + # HCubature doesn't support functions that output Unitful Quantity types + # Establish the units that are output by f + testpoint_parametriccoord = fill(FP(0.5),3) + integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) + # Create a wrapper that returns only the value component in those units + uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) + # Integrate only the unitless values + value = HCubature.hcubature(uintegrand, FP[0], FP[1]; settings.kwargs...)[1] + + # Reapply units + return value .* integrandunits +end diff --git a/src/specializations/ConeSurface.jl b/src/specializations/ConeSurface.jl new file mode 100644 index 00000000..67cbc84f --- /dev/null +++ b/src/specializations/ConeSurface.jl @@ -0,0 +1,18 @@ +################################################################################ +# Specialized Methods for ConeSurface +################################################################################ + +function integral( + f::F, + cone::Meshes.ConeSurface, + settings::I, + FP::Type{T} = Float64 +) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} + # The generic method only parameterizes the sides + sides = _integral(f, cone, settings, FP) + + # Integrate the Disk at the base + base = _integral(f, cone.base, settings, FP) + + return sides + base +end diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl new file mode 100644 index 00000000..7af8f213 --- /dev/null +++ b/src/specializations/CylinderSurface.jl @@ -0,0 +1,23 @@ +################################################################################ +# Specialized Methods for CylinderSurface +################################################################################ + +function integral( + f::F, + cyl::Meshes.CylinderSurface, + settings::I, + FP::Type{T} = Float64 +) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} + # The generic method only parameterizes the sides + sides = _integral(f, cyl, settings, FP) + + # Integrate the Disk at the top + disk_top = Meshes.Disk(cyl.top, cyl.radius) + top = _integral(f, disk_top, settings, FP) + + # Integrate the Disk at the bottom + disk_bottom = Meshes.Disk(cyl.bot, cyl.radius) + bottom = _integral(f, disk_bottom, settings, FP) + + return sides + top + bottom +end diff --git a/src/specializations/FrustumSurface.jl b/src/specializations/FrustumSurface.jl new file mode 100644 index 00000000..96bea658 --- /dev/null +++ b/src/specializations/FrustumSurface.jl @@ -0,0 +1,19 @@ +################################################################################ +# Specialized Methods for FrustumSurface +################################################################################ + +function integral( + f::F, + frust::Meshes.FrustumSurface, + settings::I, + FP::Type{T} = Float64 +) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} + # The generic method only parameterizes the sides + sides = _integral(f, frust, settings, FP) + + # Integrate the Disks at the top and bottom + top = _integral(f, frust.top, settings, FP) + bottom = _integral(f, frust.bot, settings, FP) + + return sides + top + bottom +end diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl new file mode 100644 index 00000000..0c4d431a --- /dev/null +++ b/src/specializations/Line.jl @@ -0,0 +1,70 @@ +################################################################################ +# Specialized Methods for Line +################################################################################ + +function integral( + f::F, + line::Meshes.Line, + settings::GaussLegendre, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] + xs, ws = _gausslegendre(FP, settings.n) + + # Normalize the Line s.t. line(t) is distance t from origin point + line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) + + # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) + t(x) = x / (1 - x^2) + t′(x) = (1 + x^2) / (1 - x^2)^2 + + # Integrate f along the Line + domainunits = _units(line(0)) + integrand(x) = f(line(t(x))) * t′(x) * domainunits + return sum(w .* integrand(x) for (w,x) in zip(ws, xs)) +end + +function integral( + f::F, + line::Meshes.Line, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Normalize the Line s.t. line(t) is distance t from origin point + line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) + + # Integrate f along the Line + domainunits = _units(line(0)) + integrand(t) = f(line(t)) * domainunits + return QuadGK.quadgk(integrand, FP(-Inf), FP(Inf); settings.kwargs...)[1] +end + +function integral( + f::F, + line::Meshes.Line, + settings::HAdaptiveCubature, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Normalize the Line s.t. line(t) is distance t from origin point + line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) + + # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) + t(x) = x / (1 - x^2) + t′(x) = (1 + x^2) / (1 - x^2)^2 + + # Integrate f along the Line + domainunits = _units(line(0)) + integrand(x::AbstractVector) = f(line(t(x[1]))) * t′(x[1]) * domainunits + + # HCubature doesn't support functions that output Unitful Quantity types + # Establish the units that are output by f + testpoint_parametriccoord = FP[0.5] + integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) + # Create a wrapper that returns only the value component in those units + uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) + # Integrate only the unitless values + value = HCubature.hcubature(uintegrand, FP[-1], FP[1]; settings.kwargs...)[1] + + # Reapply units + return value .* integrandunits +end diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl new file mode 100644 index 00000000..e074c5c4 --- /dev/null +++ b/src/specializations/Plane.jl @@ -0,0 +1,72 @@ +################################################################################ +# Specialized Methods for Plane +################################################################################ + +function integral( + f::F, + plane::Meshes.Plane, + settings::GaussLegendre, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² + xs, ws = _gausslegendre(FP, settings.n) + wws = Iterators.product(ws, ws) + xxs = Iterators.product(xs, xs) + + # Normalize the Plane's orthogonal vectors + plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v)) + + # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) + t(x) = x / (1 - x^2) + t′(x) = (1 + x^2) / (1 - x^2)^2 + + # Integrate f over the Plane + domainunits = _units(plane(0,0)) + integrand(((wi,wj), (xi,xj))) = wi * wj * f(plane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2 + return sum(integrand, zip(wws,xxs)) +end + +function integral( + f::F, + plane::Meshes.Plane, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Normalize the Plane's orthogonal vectors + plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v)) + + # Integrate f over the Plane + domainunits = _units(plane(0,0)) + inner∫(v) = QuadGK.quadgk(u -> f(plane(u,v)) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1] + return QuadGK.quadgk(v -> inner∫(v) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1] +end + +function integral( + f::F, + plane::Meshes.Plane, + settings::HAdaptiveCubature, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Normalize the Plane's orthogonal vectors + plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v)) + + # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) + t(x) = x / (1 - x^2) + t′(x) = (1 + x^2) / (1 - x^2)^2 + + # Integrate f over the Plane + domainunits = _units(plane(0,0)) + integrand(x::AbstractVector) = f(plane(t(x[1]), t(x[2]))) * t′(x[1]) * t′(x[2]) * domainunits^2 + + # HCubature doesn't support functions that output Unitful Quantity types + # Establish the units that are output by f + testpoint_parametriccoord = FP[0.5, 0.5] + integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) + # Create a wrapper that returns only the value component in those units + uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) + # Integrate only the unitless values + value = HCubature.hcubature(uintegrand, FP[-1,-1], FP[1,1]; settings.kwargs...)[1] + + # Reapply units + return value .* integrandunits +end diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl new file mode 100644 index 00000000..b7d3f58f --- /dev/null +++ b/src/specializations/Ray.jl @@ -0,0 +1,73 @@ +################################################################################ +# Specialized Methods for Ray +################################################################################ + +function integral( + f::F, + ray::Meshes.Ray, + settings::GaussLegendre, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] + xs, ws = _gausslegendre(FP, settings.n) + + # Normalize the Ray s.t. ray(t) is distance t from origin point + ray = Ray(ray.p, Meshes.unormalize(ray.v)) + + # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞) + t₁(x) = FP(1/2) * x + FP(1/2) + t₁′(x) = FP(1/2) + t₂(x) = x / (1 - x^2) + t₂′(x) = (1 + x^2) / (1 - x^2)^2 + t = t₂ ∘ t₁ + t′(x) = t₂′(t₁(x)) * t₁′(x) + + # Integrate f along the Ray + domainunits = _units(ray(0)) + integrand(x) = f(ray(t(x))) * t′(x) * domainunits + return sum(w .* integrand(x) for (w,x) in zip(ws, xs)) +end + +function integral( + f::F, + ray::Meshes.Ray, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Normalize the Ray s.t. ray(t) is distance t from origin point + ray = Ray(ray.p, Meshes.unormalize(ray.v)) + + # Integrate f along the Ray + domainunits = _units(ray(0)) + return QuadGK.quadgk(t -> f(ray(t)) * domainunits, FP(0), FP(Inf); settings.kwargs...)[1] +end + +function integral( + f::F, + ray::Meshes.Ray, + settings::HAdaptiveCubature, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Normalize the Ray s.t. ray(t) is distance t from origin point + ray = Ray(ray.p, Meshes.unormalize(ray.v)) + + # Domain transformation: x ∈ [0,1] ↦ t ∈ [0,∞) + t(x) = x / (1 - x^2) + t′(x) = (1 + x^2) / (1 - x^2)^2 + + # Integrate f along the Ray + domainunits = _units(ray(0)) + integrand(x::AbstractVector) = f(ray(t(x[1]))) * t′(x[1]) * domainunits + + # HCubature doesn't support functions that output Unitful Quantity types + # Establish the units that are output by f + testpoint_parametriccoord = FP[0.5] + integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) + # Create a wrapper that returns only the value component in those units + uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) + # Integrate only the unitless values + value = HCubature.hcubature(uintegrand, FP[0], FP[1]; settings.kwargs...)[1] + + # Reapply units + return value .* integrandunits +end diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl new file mode 100644 index 00000000..4876a768 --- /dev/null +++ b/src/specializations/Ring.jl @@ -0,0 +1,41 @@ +################################################################################ +# Specialized Methods for Ring +################################################################################ + +function integral( + f::F, + ring::Meshes.Ring, + settings::HAdaptiveCubature +) where {F<:Function} + # Convert the Ring into Segments, sum the integrals of those + return sum(segment -> _integral(f, segment, settings), segments(ring)) +end + +function integral( + f::F, + ring::Meshes.Ring, + settings::HAdaptiveCubature, + FP::Type{T} +) where {F<:Function, T<:AbstractFloat} + # Convert the Ring into Segments, sum the integrals of those + return sum(segment -> _integral(f, segment, settings, FP), segments(ring)) +end + +function integral( + f::F, + ring::Meshes.Ring, + settings::I +) where {F<:Function, I<:IntegrationAlgorithm} + # Convert the Ring into Segments, sum the integrals of those + return sum(segment -> integral(f, segment, settings), segments(ring)) +end + +function integral( + f::F, + ring::Meshes.Ring, + settings::I, + FP::Type{T} +) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} + # Convert the Ring into Segments, sum the integrals of those + return sum(segment -> integral(f, segment, settings, FP), segments(ring)) +end diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl new file mode 100644 index 00000000..03917618 --- /dev/null +++ b/src/specializations/Rope.jl @@ -0,0 +1,41 @@ +################################################################################ +# Specialized Methods for Rope +################################################################################ + +function integral( + f::F, + rope::Meshes.Rope, + settings::HAdaptiveCubature +) where {F<:Function} + # Convert the Rope into Segments, sum the integrals of those + return sum(segment -> _integral(f, segment, settings), segments(rope)) +end + +function integral( + f::F, + rope::Meshes.Rope, + settings::HAdaptiveCubature, + FP::Type{T} +) where {F<:Function, T<:AbstractFloat} + # Convert the Rope into Segments, sum the integrals of those + return sum(segment -> _integral(f, segment, settings, FP), segments(rope)) +end + +function integral( + f::F, + rope::Meshes.Rope, + settings::I +) where {F<:Function, I<:IntegrationAlgorithm} + # Convert the Rope into Segments, sum the integrals of those + return sum(segment -> integral(f, segment, settings), segments(rope)) +end + +function integral( + f::F, + rope::Meshes.Rope, + settings::I, + FP::Type{T} +) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat} + # Convert the Rope into Segments, sum the integrals of those + return sum(segment -> integral(f, segment, settings, FP), segments(rope)) +end diff --git a/src/integral_volume.jl b/src/specializations/Tetrahedron.jl similarity index 53% rename from src/integral_volume.jl rename to src/specializations/Tetrahedron.jl index 60943a86..51beec43 100644 --- a/src/integral_volume.jl +++ b/src/specializations/Tetrahedron.jl @@ -1,40 +1,3 @@ -################################################################################ -# Generalized 3D Methods -################################################################################ - -function _integral_3d( - f, - geometry3d, - settings::GaussLegendre, - FP::Type{T} = Float64 -) where {T<:AbstractFloat} - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 - xs, ws = _gausslegendre(FP, settings.n) - wws = Iterators.product(ws, ws, ws) - xxs = Iterators.product(xs, xs, xs) - - # Domain transformation: x [-1,1] ↦ s,t,u [0,1] - t(x) = FP(1/2) * x + FP(1/2) - - function integrand(((wi,wj,wk), (xi,xj,xk))) - ts = t.([xi, xj, xk]) - wi * wj * wk * f(geometry3d(ts...)) * differential(geometry3d, ts) - end - - return FP(1/8) .* sum(integrand, zip(wws,xxs)) -end - -# Integrating volumes with GaussKronrod not supported by default -function _integral_3d( - f, - geometry, - settings::GaussKronrod, - FP::Type{T} = Float64 -) where {T<:AbstractFloat} - error("Integrating this volume type with GaussKronrod not supported.") -end - - ################################################################################ # Specialized Methods for Tetrahedron ################################################################################ diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl new file mode 100644 index 00000000..53982dc9 --- /dev/null +++ b/src/specializations/Triangle.jl @@ -0,0 +1,99 @@ +################################################################################ +# Specialized Methods for Triangle +################################################################################ + +""" + integral(f, triangle::Meshes.Triangle, ::GaussLegendre) + +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 integral( + f::F, + triangle::Meshes.Ngon{3}, + settings::GaussLegendre, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 + xs, ws = _gausslegendre(FP, settings.n) + wws = Iterators.product(ws, ws) + xxs = Iterators.product(xs, xs) + + # Domain transformations: + # xᵢ [-1,1] ↦ R [0,1] + # xⱼ [-1,1] ↦ φ [0,π/2] + uR(xᵢ) = T(1/2) * (xᵢ + 1) + uφ(xⱼ) = T(π/4) * (xⱼ + 1) + + # Integrate the Barycentric triangle by transforming it into polar coordinates + # with a modified radius + # R = r ( sinφ + cosφ ) + # s.t. integration bounds become rectangular + # R ∈ [0, 1] and φ ∈ [0, π/2] + function integrand(((wᵢ,wⱼ), (xᵢ,xⱼ))) + R = uR(xᵢ) + φ = uφ(xⱼ) + a,b = sincos(φ) + u = R * (1 - a / (a + b)) + v = R * (1 - b / (a + b)) + return wᵢ * wⱼ * f(triangle(u, v)) * R / (a + b)^2 + end + + # Calculate 2D Gauss-Legendre integral over modified-polar-Barycentric coordinates + # Apply a linear domain-correction factor + return FP(π/4) * area(triangle) .* sum(integrand, zip(wws,xxs)) +end + +""" + integral(f, triangle::Meshes.Triangle, ::GaussKronrod) + +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 integral( + f::F, + triangle::Meshes.Ngon{3}, + settings::GaussKronrod, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0) + # i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du + inner∫(u) = QuadGK.quadgk(v -> f(triangle(u,v)), FP(0), FP(1-u); settings.kwargs...)[1] + outer∫ = QuadGK.quadgk(inner∫, FP(0), FP(1); settings.kwargs...)[1] + + # Apply a linear domain-correction factor 0.5 ↦ area(triangle) + return 2 * area(triangle) .* outer∫ +end + +""" + integral(f, triangle::Meshes.Triangle, ::GaussKronrod) + +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 integral( + f::F, + triangle::Meshes.Ngon{3}, + settings::HAdaptiveCubature, + FP::Type{T} = Float64 +) where {F<:Function, T<:AbstractFloat} + # Integrate the Barycentric triangle by transforming it into polar coordinates + # with a modified radius + # R = r ( sinφ + cosφ ) + # s.t. integration bounds become rectangular + # R ∈ [0, 1] and φ ∈ [0, π/2] + function integrand(Rφ) + R,φ = Rφ + a,b = sincos(φ) + u = R * (1 - a / (a + b)) + v = R * (1 - b / (a + b)) + return f(triangle(u, v)) * R / (a + b)^2 + end + intval = HCubature.hcubature(integrand, FP[0, 0], FP[1, π/2], settings.kwargs...)[1] + + # Apply a linear domain-correction factor 0.5 ↦ area(triangle) + return 2 * area(triangle) .* intval +end diff --git a/src/utils.jl b/src/utils.jl index 45782e82..89259b35 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -10,9 +10,11 @@ central-finite-difference approximation with step size `ε`. """ function jacobian( geometry, - ts::AbstractVector{T}; + ts; ε=1e-6 -) where {T<:AbstractFloat} +) + T = eltype(ts) + # Get the partial derivative along the n'th axis via finite difference approximation # where ts is the current parametric position (εv is a reusable buffer) function ∂ₙr!(εv, ts, n) @@ -27,33 +29,30 @@ function jacobian( # Central finite difference function ∂ₙr_central!(εv, ts, n) - εv .= T(0) εv[n] = ε - a = ts - εv - b = ts + εv + a = ts .- εv + b = ts .+ εv return (geometry(b...) - geometry(a...)) / 2ε end # Left finite difference function ∂ₙr_left!(εv, ts, n) - εv .= T(0) εv[n] = ε - a = ts - εv + a = ts .- εv b = ts return (geometry(b...) - geometry(a...)) / ε end # Right finite difference function ∂ₙr_right!(εv, ts, n) - εv .= T(0) εv[n] = ε a = ts - b = ts + εv + b = ts .+ εv return (geometry(b...) - geometry(a...)) / ε end # Allocate a re-usable ε vector - εv = similar(ts) + εv = zeros(T, length(ts)) ∂ₙr(n) = ∂ₙr!(εv,ts,n) return map(∂ₙr, 1:length(ts)) @@ -65,7 +64,10 @@ end Calculate the differential element (length, area, volume, etc) of the parametric function for `geometry` at arguments `ts`. """ -function differential(geometry, ts::AbstractVector) +function differential( + geometry, + ts +) J = jacobian(geometry, ts) # TODO generalize this with geometric algebra, e.g.: norm(foldl(∧, J)) diff --git a/test/Project.toml b/test/Project.toml index 97be0fa2..8ba08814 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,5 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -MeshIntegrals = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47" Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -10,7 +9,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] Aqua = "0.8" -MeshIntegrals = "0.13" Meshes = "0.50, 0.51" QuadGK = "2" TestItemRunner = "1" diff --git a/test/combinations.jl b/test/combinations.jl index 84d122e4..0707bac6 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -50,14 +50,14 @@ end # 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, 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 - @test integral(fv, cone, GaussKronrod()) ≈ vsol + @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 @@ -94,14 +94,14 @@ end area_walls = area_walls_projected - area_walls_missing area_walls + _area_base(r_top) + _area_base(r_bot) end - @test integral(f, frustum, GaussLegendre(100)) ≈ sol - @test integral(f, frustum, GaussKronrod()) ≈ sol + @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 - @test integral(fv, frustum, GaussKronrod()) ≈ vsol + @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 end