From ce097f66c43a2f39c96e832634d0c50f83b06442 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:35:30 -0400 Subject: [PATCH 01/89] Add JacobianMethod types --- src/differentiation.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index acb580f1..886e4120 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -1,5 +1,21 @@ ################################################################################ -# Jacobian and Differential Elements +# JacobianMethods +################################################################################ + +abstract type JacobianMethod end + +struct FiniteDifference{T<:AbstractFloat} <: JacobianMethod + ε::T +end + +# struct EnzymeAD <: JacobianMethod +# end + +# struct ZygoteAD <: JacobianMethod +# end + +################################################################################ +# Jacobian ################################################################################ """ From 3decf9f4354a303707510265e3dbddd536bbce1f Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:43:10 -0400 Subject: [PATCH 02/89] Change JacobianMethod to DifferentiationMethod --- src/differentiation.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 886e4120..4bf95a3b 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -2,16 +2,22 @@ # JacobianMethods ################################################################################ -abstract type JacobianMethod end +abstract type DifferentiationMethod end -struct FiniteDifference{T<:AbstractFloat} <: JacobianMethod +""" + FiniteDifference(ε) + +Use a finite-difference approximation method to calculate derivatives with a +step size of `ε`. +""" +struct FiniteDifference{T<:AbstractFloat} <: DifferentiationMethod ε::T end -# struct EnzymeAD <: JacobianMethod +# struct EnzymeAD <: DifferentiationMethod # end -# struct ZygoteAD <: JacobianMethod +# struct ZygoteAD <: DifferentiationMethod # end ################################################################################ From bd54b312745c8f64c7352bd1930c99f5a8a86cb1 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:44:18 -0400 Subject: [PATCH 03/89] Add a third argument to jacobian for DiffMethod for dispatch --- src/differentiation.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 4bf95a3b..d80e22c3 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -25,19 +25,22 @@ end ################################################################################ """ - jacobian(geometry, ts; ε=1e-6) + jacobian(geometry, ts, method) -Calculate the Jacobian of a geometry at some parametric point `ts` using a simple -finite-difference approximation with step size `ε`. +Calculate the Jacobian of a geometry's parametric function at some point `ts` +using a particular differentiation method. # Arguments - `geometry`: some `Meshes.Geometry` of N parametric dimensions - `ts`: a parametric point specified as a vector or tuple of length N -- `ε`: step size to use for the finite-difference approximation +- `method`: the desired `DifferentiationMethod` """ +function jacobian end + function jacobian( geometry::Geometry, - ts::V; + ts::V, + fd::FiniteDifference; ε = 1e-6 ) where {V <: Union{AbstractVector, Tuple}} Dim = Meshes.paramdim(geometry) @@ -46,7 +49,7 @@ function jacobian( end T = eltype(ts) - ε = T(ε) + ε = T(fd.ε) # Get the partial derivative along the n'th axis via finite difference # approximation, where ts is the current parametric position From 4cb6b0224c61e3e480ab385ff8747ecca6c8d385 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:45:32 -0400 Subject: [PATCH 04/89] Remove dead line --- src/differentiation.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index d80e22c3..df66677b 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -40,8 +40,7 @@ function jacobian end function jacobian( geometry::Geometry, ts::V, - fd::FiniteDifference; - ε = 1e-6 + fd::FiniteDifference ) where {V <: Union{AbstractVector, Tuple}} Dim = Meshes.paramdim(geometry) if Dim != length(ts) From 34fbca2e092b0faa854354f49a68e1ab4387095c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:45:51 -0400 Subject: [PATCH 05/89] Remove obsolete analytic jacobian for BezierCurve --- src/differentiation.jl | 30 ------------------------------ 1 file changed, 30 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index df66677b..9071815d 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -72,36 +72,6 @@ function jacobian( return ntuple(n -> ∂ₙr(ts, n), Dim) end -function jacobian( - bz::Meshes.BezierCurve, - ts::V -) where {V <: Union{AbstractVector, Tuple}} - t = only(ts) - # Parameter t restricted to domain [0,1] by definition - if t < 0 || t > 1 - throw(DomainError(t, "b(t) is not defined for t outside [0, 1].")) - end - - # Aliases - P = bz.controls - N = Meshes.degree(bz) - - # Ensure that this implementation is tractible: limited by ability to calculate - # binomial(N, N/2) without overflow. It's possible to extend this range by - # converting N to a BigInt, but this results in always returning BigFloat's. - N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.") - - # Generator for Bernstein polynomial functions - B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i) - - # Derivative = N Σ_{i=0}^{N-1} sigma(i) - # P indices adjusted for Julia 1-based array indexing - sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1]) - derivative = N .* sum(sigma, 0:(N - 1)) - - return (derivative,) -end - ################################################################################ # Differential Elements ################################################################################ From 7e5a405988afe6effddec9789e1e994cf34112d8 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:50:20 -0400 Subject: [PATCH 06/89] Add default step size if unspecified --- src/differentiation.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/differentiation.jl b/src/differentiation.jl index 9071815d..40cecb9a 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -14,6 +14,9 @@ struct FiniteDifference{T<:AbstractFloat} <: DifferentiationMethod ε::T end +# If ε not specified, default to 1e-6 +FiniteDifference() = new(1e-6) + # struct EnzymeAD <: DifferentiationMethod # end From f17f326237473664deed874256e03bff2c06e960 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 22:51:01 -0400 Subject: [PATCH 07/89] Convert differential to three-argument with default FiniteDiff method --- src/differentiation.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 40cecb9a..25ee153c 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -80,17 +80,18 @@ end ################################################################################ """ - differential(geometry, ts) + differential(geometry, ts, method=FiniteDifference()) Calculate the differential element (length, area, volume, etc) of the parametric function for `geometry` at arguments `ts`. """ function differential( geometry::Geometry, - ts::V + ts::V, + method::DifferentiationMethod = FiniteDifference() ) where {V <: Union{AbstractVector, Tuple}} # Calculate the Jacobian, convert Vec -> KVector - J = jacobian(geometry, ts) + J = jacobian(geometry, ts, method) J_kvecs = Iterators.map(_kvector, J) # Extract units from Geometry type From ac23e4335e9f4d7feeca9adef0c5da1dd55d972e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 23:20:41 -0400 Subject: [PATCH 08/89] Bigfix --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 25ee153c..617724b4 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -15,7 +15,7 @@ struct FiniteDifference{T<:AbstractFloat} <: DifferentiationMethod end # If ε not specified, default to 1e-6 -FiniteDifference() = new(1e-6) +FiniteDifference() = FiniteDifference(1e-6) # struct EnzymeAD <: DifferentiationMethod # end From ceca1b5a1ef3c5ae521bd56520d943941469a8c8 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 25 Oct 2024 23:25:37 -0400 Subject: [PATCH 09/89] Add two-arg Jacobian pointing to FD --- src/differentiation.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/differentiation.jl b/src/differentiation.jl index 617724b4..f8e07b40 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -40,6 +40,13 @@ using a particular differentiation method. """ function jacobian end +function jacobian( + geometry::Geometry, + ts::V +) where {V <: Union{AbstractVector, Tuple}} + return jacobian(geometry, ts, FiniteDifference()) +end + function jacobian( geometry::Geometry, ts::V, From 2e7b23f18ceffc9e68f4e78a4f09b3094cdee71c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 01:14:22 -0400 Subject: [PATCH 10/89] Apply formatting suggestion Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index f8e07b40..008b33b2 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -10,7 +10,7 @@ abstract type DifferentiationMethod end Use a finite-difference approximation method to calculate derivatives with a step size of `ε`. """ -struct FiniteDifference{T<:AbstractFloat} <: DifferentiationMethod +struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod ε::T end From bfab8836216fc338e76d61a10eb0927161371c59 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 07:33:08 -0400 Subject: [PATCH 11/89] Apply suggestion Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 008b33b2..c70a931f 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -28,7 +28,7 @@ FiniteDifference() = FiniteDifference(1e-6) ################################################################################ """ - jacobian(geometry, ts, method) + jacobian(geometry, ts, method=FiniteDifference()) Calculate the Jacobian of a geometry's parametric function at some point `ts` using a particular differentiation method. From 62d7a531e4fa260d56d76a726233ba1cb8c15e4c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 08:43:36 -0400 Subject: [PATCH 12/89] Add a dt argument --- src/integral.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index d4925db7..ba7b30a7 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -59,8 +59,9 @@ function _integral( f, geometry, rule::GaussLegendre; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {T <: AbstractFloat} +) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) # Get Gauss-Legendre nodes and weights for a region [-1,1]^N @@ -73,7 +74,7 @@ function _integral( function integrand((weights, nodes)) ts = t.(nodes) - prod(weights) * f(geometry(ts...)) * differential(geometry, ts) + prod(weights) * f(geometry(ts...)) * differential(geometry, ts, dt) end return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes)) @@ -84,11 +85,12 @@ function _integral( f, geometry, rule::HAdaptiveCubature; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {T <: AbstractFloat} N = Meshes.paramdim(geometry) - integrand(ts) = f(geometry(ts...)) * differential(geometry, ts) + integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, dt) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f @@ -111,9 +113,10 @@ function _integral_gk_1d( f, geometry, rule::GaussKronrod; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {T <: AbstractFloat} - integrand(t) = f(geometry(t)) * differential(geometry, (t,)) + integrand(t) = f(geometry(t)) * differential(geometry, (t,), dt) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end @@ -121,9 +124,10 @@ function _integral_gk_2d( f, geometry2d, rule::GaussKronrod; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {T <: AbstractFloat} - integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v)) + integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), dt) ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1] end From 3114e549798fecfec0da9661e12f0ed489abfd24 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 08:45:58 -0400 Subject: [PATCH 13/89] Update docstring API --- src/integral.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/integral.jl b/src/integral.jl index ba7b30a7..65ba1420 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - integral(f, geometry[, rule]; FP=Float64) + integral(f, geometry[, rule]; dt=FiniteDifference(), FP=Float64) Numerically integrate a given function `f(::Point)` over the domain defined by a `geometry` using a particular numerical integration `rule` with floating point @@ -14,6 +14,7 @@ precision of type `FP`. - `geometry`: some `Meshes.Geometry` that defines the integration domain - `rule`: optionally, the `IntegrationRule` used for integration (by default `GaussKronrod()` in 1D and `HAdaptiveCubature()` else) - `FP`: optionally, the floating point precision desired (`Float64` by default) +- `dt`: optionally, use a particular `DifferentiationMethod` for calculating differential elements Note that reducing `FP` below `Float64` will incur some loss of precision. By contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision From 3f0df9aa90553ce436bc82bc949a41e2850ee0ea Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 08:50:09 -0400 Subject: [PATCH 14/89] Add dt kwargs --- src/specializations/BezierCurve.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index f0088e5f..8e167f4d 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -11,7 +11,7 @@ """ integral(f, curve::BezierCurve, ::GaussLegendre; - FP=Float64, alg=Meshes.Horner()) + dt=FiniteDifference(), FP=Float64, 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 @@ -24,16 +24,17 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussLegendre; + dt::DM = FiniteDifference(), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] xs, ws = _gausslegendre(FP, rule.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) - integrand(x) = f(point(x)) * differential(curve, (t(x),)) + integrand(x) = f(point(x)) * differential(curve, (t(x),), dt) # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length] return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs)) @@ -41,7 +42,7 @@ end """ integral(f, curve::BezierCurve, ::GaussKronrod; - FP=Float64, alg=Meshes.Horner()) + dt=DifferentiationMethod(), FP=Float64, 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 @@ -54,17 +55,18 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussKronrod; + dt::DM = FiniteDifference(), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} point(t) = curve(t, alg) - integrand(t) = f(point(t)) * differential(curve, (t,)) + integrand(t) = f(point(t)) * differential(curve, (t,), dt) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end """ integral(f, curve::BezierCurve, ::HAdaptiveCubature; - FP=Float64, alg=Meshes.Horner()) + dt=FiniteDifference(), FP=Float64, 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 @@ -77,11 +79,12 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::HAdaptiveCubature; + dt::DM = FiniteDifference(), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} point(t) = curve(t, alg) - integrand(ts) = f(point(only(ts))) * differential(curve, ts) + integrand(ts) = f(point(only(ts))) * differential(curve, ts, dt) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f From 8d7753dfbfbdb7e376fae13f0216c6f339687280 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 08:51:14 -0400 Subject: [PATCH 15/89] Add missing params --- src/integral.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 65ba1420..4a45db88 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -88,7 +88,7 @@ function _integral( rule::HAdaptiveCubature; dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {T <: AbstractFloat} +) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, dt) @@ -116,7 +116,7 @@ function _integral_gk_1d( rule::GaussKronrod; dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {T <: AbstractFloat} +) where {DM <: DifferentiationMethod, T <: AbstractFloat} integrand(t) = f(geometry(t)) * differential(geometry, (t,), dt) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end @@ -127,7 +127,7 @@ function _integral_gk_2d( rule::GaussKronrod; dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {T <: AbstractFloat} +) where {DM <: DifferentiationMethod, T <: AbstractFloat} integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), dt) ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1] From 855f3a9ef3a57bc800aac26150dfe402b5bf6560 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 08:57:21 -0400 Subject: [PATCH 16/89] Add dt kwargs (with TODOs) --- src/specializations/Line.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index fbb398c9..04fbf2c2 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -11,8 +11,9 @@ function integral( f::F, line::Meshes.Line, rule::GaussLegendre; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] xs, ws = _gausslegendre(FP, rule.n) @@ -25,7 +26,7 @@ function integral( # Integrate f along the Line domainunits = _units(line(0)) - integrand(x) = f(line(t(x))) * t′(x) * domainunits + integrand(x) = f(line(t(x))) * t′(x) * domainunits # TODO differential return sum(w .* integrand(x) for (w, x) in zip(ws, xs)) end @@ -33,14 +34,15 @@ function integral( f::F, line::Meshes.Line, rule::GaussKronrod; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Line s.t. line(t) is distance t from origin point line = Meshes.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 + integrand(t) = f(line(t)) * domainunits # TODO differential return QuadGK.quadgk(integrand, FP(-Inf), FP(Inf); rule.kwargs...)[1] end @@ -48,8 +50,9 @@ function integral( f::F, line::Meshes.Line, rule::HAdaptiveCubature; + dt::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Line s.t. line(t) is distance t from origin point line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) @@ -59,7 +62,7 @@ function integral( # Integrate f along the Line domainunits = _units(line(0)) - integrand(x::AbstractVector) = f(line(t(x[1]))) * t′(x[1]) * domainunits + integrand(x::AbstractVector) = f(line(t(x[1]))) * t′(x[1]) * domainunits # TODO differential # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f From d832b0589911f483b4d2741ef89878e96bbb058a Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 09:16:10 -0400 Subject: [PATCH 17/89] Consolidate kwarg naming --- src/differentiation.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index c70a931f..a7f4eff5 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -28,15 +28,15 @@ FiniteDifference() = FiniteDifference(1e-6) ################################################################################ """ - jacobian(geometry, ts, method=FiniteDifference()) + jacobian(geometry, ts, dt=FiniteDifference()) Calculate the Jacobian of a geometry's parametric function at some point `ts` -using a particular differentiation method. +using a particular differentiation method `dt`. # Arguments - `geometry`: some `Meshes.Geometry` of N parametric dimensions - `ts`: a parametric point specified as a vector or tuple of length N -- `method`: the desired `DifferentiationMethod` +- `dt`: the desired `DifferentiationMethod` """ function jacobian end @@ -50,7 +50,7 @@ end function jacobian( geometry::Geometry, ts::V, - fd::FiniteDifference + dt::FiniteDifference ) where {V <: Union{AbstractVector, Tuple}} Dim = Meshes.paramdim(geometry) if Dim != length(ts) @@ -58,7 +58,7 @@ function jacobian( end T = eltype(ts) - ε = T(fd.ε) + ε = T(dt.ε) # Get the partial derivative along the n'th axis via finite difference # approximation, where ts is the current parametric position @@ -87,7 +87,7 @@ end ################################################################################ """ - differential(geometry, ts, method=FiniteDifference()) + differential(geometry, ts, dt=FiniteDifference()) Calculate the differential element (length, area, volume, etc) of the parametric function for `geometry` at arguments `ts`. @@ -95,10 +95,10 @@ function for `geometry` at arguments `ts`. function differential( geometry::Geometry, ts::V, - method::DifferentiationMethod = FiniteDifference() + dt::DifferentiationMethod = FiniteDifference() ) where {V <: Union{AbstractVector, Tuple}} # Calculate the Jacobian, convert Vec -> KVector - J = jacobian(geometry, ts, method) + J = jacobian(geometry, ts, dt) J_kvecs = Iterators.map(_kvector, J) # Extract units from Geometry type From 6935da6427349cf24149e44b3e12b7c979af12e7 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 09:16:29 -0400 Subject: [PATCH 18/89] Update src/differentiation.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index a7f4eff5..eed5d407 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -5,7 +5,7 @@ abstract type DifferentiationMethod end """ - FiniteDifference(ε) + FiniteDifference(ε=1e-6) Use a finite-difference approximation method to calculate derivatives with a step size of `ε`. From bfcd204ef7878cd7f8063a7b120a45cb4a57e367 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 09:57:45 -0400 Subject: [PATCH 19/89] Add exports --- src/MeshIntegrals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 19228769..cbd6c158 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -12,7 +12,7 @@ import Unitful include("utils.jl") include("differentiation.jl") -export jacobian +export DifferentiationMethod, FiniteDifference, jacobian include("integration_rules.jl") export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature From 80b05f99aa461b9e736eb94e3740d1bb9330da22 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 10:49:40 -0400 Subject: [PATCH 20/89] Tweak suggested naming [skip CI] --- src/differentiation.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index eed5d407..ccce7c73 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -17,11 +17,9 @@ end # If ε not specified, default to 1e-6 FiniteDifference() = FiniteDifference(1e-6) -# struct EnzymeAD <: DifferentiationMethod -# end +# struct AutoEnzyme <: DifferentiationMethod end -# struct ZygoteAD <: DifferentiationMethod -# end +# struct AutoZygote <: DifferentiationMethod end ################################################################################ # Jacobian From 323c5d0325b19d820c5e025f6228197bcce8755e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 11:08:43 -0400 Subject: [PATCH 21/89] Update Line methods to use an internal differential method --- src/specializations/Line.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 04fbf2c2..87a06c76 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -25,8 +25,8 @@ function integral( 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 # TODO differential + differential(line, x) = t′(x) * _units(line(0)) + integrand(x) = f(line(t(x))) * differential(line, x) return sum(w .* integrand(x) for (w, x) in zip(ws, xs)) end @@ -61,8 +61,8 @@ function integral( 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 # TODO differential + differential(line, x) = t′(x) * _units(line(0)) + integrand(x::AbstractVector) = f(line(t(x[1]))) * differential(line, x[1]) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f @@ -71,7 +71,7 @@ function integral( # 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]; rule.kwargs...)[1] + value = HCubature.hcubature(uintegrand, (-one(FP),), (one(FP),); rule.kwargs...)[1] # Reapply units return value .* integrandunits From 25ce37d3f5e60da6130a32f62d4277f57f4e9677 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 13:44:24 -0400 Subject: [PATCH 22/89] Bump version to v0.16.0-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7e9cfe3b..5a1d4711 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MeshIntegrals" uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47" -version = "0.15.2" +version = "0.16.0-DEV" [deps] CliffordNumbers = "3998ac73-6bd4-4031-8035-f167dd3ed523" From 78c35629ec1284bafe59da167e35b968255d9bb6 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 13:56:45 -0400 Subject: [PATCH 23/89] Add block for DifferentiationMethods --- docs/src/api.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/src/api.md b/docs/src/api.md index 0d231931..40cc8799 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -24,6 +24,11 @@ MeshIntegrals.HAdaptiveCubature ## Derivatives +```@docs +MeshIntegrals.DifferentiationMethod +MeshIntegrals.FiniteDifference +``` + ```@docs MeshIntegrals.jacobian ``` From 855d6c6739adaacf7cbe526a9d8202a37f6aeb94 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 13:57:31 -0400 Subject: [PATCH 24/89] Drop abstract type --- docs/src/api.md | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/src/api.md b/docs/src/api.md index 40cc8799..27888fc8 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -25,7 +25,6 @@ MeshIntegrals.HAdaptiveCubature ## Derivatives ```@docs -MeshIntegrals.DifferentiationMethod MeshIntegrals.FiniteDifference ``` From 109193c1ecf12ef3a05d864e548f574080c79b98 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 15:16:30 -0400 Subject: [PATCH 25/89] Add diff_method arguments --- src/specializations/Plane.jl | 9 ++++++--- src/specializations/Ray.jl | 9 ++++++--- src/specializations/Tetrahedron.jl | 9 ++++++--- src/specializations/Triangle.jl | 9 ++++++--- 4 files changed, 24 insertions(+), 12 deletions(-) diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 70243e8a..79039868 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -12,8 +12,9 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussLegendre; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² xs, ws = _gausslegendre(FP, rule.n) wws = Iterators.product(ws, ws) @@ -40,8 +41,9 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussKronrod; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Plane's orthogonal vectors uu = Meshes.unormalize(plane.u) uv = Meshes.unormalize(plane.v) @@ -58,8 +60,9 @@ function integral( f::F, plane::Meshes.Plane, rule::HAdaptiveCubature; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Plane's orthogonal vectors uu = Meshes.unormalize(plane.u) uv = Meshes.unormalize(plane.v) diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index fce24530..e17f64a6 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -12,8 +12,9 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussLegendre; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] xs, ws = _gausslegendre(FP, rule.n) @@ -38,8 +39,9 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussKronrod; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Ray s.t. ray(t) is distance t from origin point ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) @@ -53,8 +55,9 @@ function integral( f::F, ray::Meshes.Ray, rule::HAdaptiveCubature; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Ray s.t. ray(t) is distance t from origin point ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index e4aac959..4c9c52d3 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -13,8 +13,9 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussLegendre; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _error_unsupported_combination("Tetrahedron", "GaussLegendre") end @@ -22,8 +23,9 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussKronrod; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} nil = zero(FP) ∫uvw(u, v, w) = f(tetrahedron(u, v, w)) ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), nil, FP(1 - v - w); rule.kwargs...)[1] @@ -38,7 +40,8 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::HAdaptiveCubature; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature") end diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index d0b418bb..acca16b2 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -21,8 +21,9 @@ function integral( f::F, triangle::Meshes.Ngon{3}, rule::GaussLegendre; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 xs, ws = _gausslegendre(FP, rule.n) wws = Iterators.product(ws, ws) @@ -63,8 +64,9 @@ function integral( f::F, triangle::Meshes.Ngon{3}, rule::GaussKronrod; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, 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 ∫u(u) = QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1] @@ -85,8 +87,9 @@ function integral( f::F, triangle::Meshes.Ngon{3}, rule::HAdaptiveCubature; + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Integrate the Barycentric triangle by transforming it into polar coordinates # with a modified radius # R = r ( sinφ + cosφ ) From 73219613cb37f30aca6311ab0bd34c45811d9468 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 26 Oct 2024 15:16:51 -0400 Subject: [PATCH 26/89] Change dt to diff_method --- src/differentiation.jl | 16 ++++++++-------- src/integral.jl | 24 +++++++++++++----------- src/specializations/BezierCurve.jl | 18 +++++++++--------- src/specializations/Line.jl | 6 +++--- 4 files changed, 33 insertions(+), 31 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index ccce7c73..c8554188 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -26,15 +26,15 @@ FiniteDifference() = FiniteDifference(1e-6) ################################################################################ """ - jacobian(geometry, ts, dt=FiniteDifference()) + jacobian(geometry, ts, diff_method=FiniteDifference()) Calculate the Jacobian of a geometry's parametric function at some point `ts` -using a particular differentiation method `dt`. +using a particular differentiation method `diff_method`. # Arguments - `geometry`: some `Meshes.Geometry` of N parametric dimensions - `ts`: a parametric point specified as a vector or tuple of length N -- `dt`: the desired `DifferentiationMethod` +- `diff_method`: the desired `DifferentiationMethod` """ function jacobian end @@ -48,7 +48,7 @@ end function jacobian( geometry::Geometry, ts::V, - dt::FiniteDifference + diff_method::FiniteDifference ) where {V <: Union{AbstractVector, Tuple}} Dim = Meshes.paramdim(geometry) if Dim != length(ts) @@ -56,7 +56,7 @@ function jacobian( end T = eltype(ts) - ε = T(dt.ε) + ε = T(diff_method.ε) # Get the partial derivative along the n'th axis via finite difference # approximation, where ts is the current parametric position @@ -85,7 +85,7 @@ end ################################################################################ """ - differential(geometry, ts, dt=FiniteDifference()) + differential(geometry, ts, diff_method=FiniteDifference()) Calculate the differential element (length, area, volume, etc) of the parametric function for `geometry` at arguments `ts`. @@ -93,10 +93,10 @@ function for `geometry` at arguments `ts`. function differential( geometry::Geometry, ts::V, - dt::DifferentiationMethod = FiniteDifference() + diff_method::DifferentiationMethod = FiniteDifference() ) where {V <: Union{AbstractVector, Tuple}} # Calculate the Jacobian, convert Vec -> KVector - J = jacobian(geometry, ts, dt) + J = jacobian(geometry, ts, diff_method) J_kvecs = Iterators.map(_kvector, J) # Extract units from Geometry type diff --git a/src/integral.jl b/src/integral.jl index 4a45db88..fa5032f9 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - integral(f, geometry[, rule]; dt=FiniteDifference(), FP=Float64) + integral(f, geometry[, rule]; diff_method=FiniteDifference(), FP=Float64) Numerically integrate a given function `f(::Point)` over the domain defined by a `geometry` using a particular numerical integration `rule` with floating point @@ -12,9 +12,11 @@ precision of type `FP`. # Arguments - `f`: an integrand function with a method `f(::Meshes.Point)` - `geometry`: some `Meshes.Geometry` that defines the integration domain -- `rule`: optionally, the `IntegrationRule` used for integration (by default `GaussKronrod()` in 1D and `HAdaptiveCubature()` else) +- `rule`: optionally, the `IntegrationRule` used for integration (by default +`GaussKronrod()` in 1D and `HAdaptiveCubature()` else) - `FP`: optionally, the floating point precision desired (`Float64` by default) -- `dt`: optionally, use a particular `DifferentiationMethod` for calculating differential elements +- `diff_method`: optionally, use a particular `DifferentiationMethod` for +calculating differential elements Note that reducing `FP` below `Float64` will incur some loss of precision. By contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision @@ -60,7 +62,7 @@ function _integral( f, geometry, rule::GaussLegendre; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) @@ -75,7 +77,7 @@ function _integral( function integrand((weights, nodes)) ts = t.(nodes) - prod(weights) * f(geometry(ts...)) * differential(geometry, ts, dt) + prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method) end return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes)) @@ -86,12 +88,12 @@ function _integral( f, geometry, rule::HAdaptiveCubature; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) - integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, dt) + integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, diff_method) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f @@ -114,10 +116,10 @@ function _integral_gk_1d( f, geometry, rule::GaussKronrod; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - integrand(t) = f(geometry(t)) * differential(geometry, (t,), dt) + integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end @@ -125,10 +127,10 @@ function _integral_gk_2d( f, geometry2d, rule::GaussKronrod; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), dt) + integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method) ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1] end diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 8e167f4d..5be80a63 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -11,7 +11,7 @@ """ integral(f, curve::BezierCurve, ::GaussLegendre; - dt=FiniteDifference(), FP=Float64, alg=Meshes.Horner()) + diff_method=FiniteDifference(), FP=Float64, 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 @@ -24,7 +24,7 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussLegendre; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} @@ -34,7 +34,7 @@ function integral( # Change of variables: x [-1,1] ↦ t [0,1] t(x) = FP(1 // 2) * x + FP(1 // 2) point(x) = curve(t(x), alg) - integrand(x) = f(point(x)) * differential(curve, (t(x),), dt) + integrand(x) = f(point(x)) * differential(curve, (t(x),), diff_method) # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length] return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs)) @@ -42,7 +42,7 @@ end """ integral(f, curve::BezierCurve, ::GaussKronrod; - dt=DifferentiationMethod(), FP=Float64, alg=Meshes.Horner()) + diff_method=DifferentiationMethod(), FP=Float64, 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 @@ -55,18 +55,18 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussKronrod; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} point(t) = curve(t, alg) - integrand(t) = f(point(t)) * differential(curve, (t,), dt) + integrand(t) = f(point(t)) * differential(curve, (t,), diff_method) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end """ integral(f, curve::BezierCurve, ::HAdaptiveCubature; - dt=FiniteDifference(), FP=Float64, alg=Meshes.Horner()) + diff_method=FiniteDifference(), FP=Float64, 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 @@ -79,12 +79,12 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::HAdaptiveCubature; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} point(t) = curve(t, alg) - integrand(ts) = f(point(only(ts))) * differential(curve, ts, dt) + integrand(ts) = f(point(only(ts))) * differential(curve, ts, diff_method) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 87a06c76..25617474 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -11,7 +11,7 @@ function integral( f::F, line::Meshes.Line, rule::GaussLegendre; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] @@ -34,7 +34,7 @@ function integral( f::F, line::Meshes.Line, rule::GaussKronrod; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Line s.t. line(t) is distance t from origin point @@ -50,7 +50,7 @@ function integral( f::F, line::Meshes.Line, rule::HAdaptiveCubature; - dt::DM = FiniteDifference(), + diff_method::DM = FiniteDifference(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Line s.t. line(t) is distance t from origin point From 8bd3a3d1184070761959036a28a9e00bfe12a7d7 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 27 Oct 2024 12:46:15 -0400 Subject: [PATCH 27/89] Prepare to revert jacobian of BezierCurve (pending new benchmarks) --- src/specializations/BezierCurve.jl | 40 ++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 08f653e1..d393797a 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -9,6 +9,9 @@ # keyword argument and pass through the specified algorithm choice. ################################################################################ +################################################################################ +# integral +################################################################################ """ integral(f, curve::BezierCurve, rule = GaussKronrod(); diff_method=FiniteDifference(), FP=Float64, alg=Meshes.Horner()) @@ -76,3 +79,40 @@ function integral( # Reapply units return value .* integrandunits end + +################################################################################ +# jacobian +################################################################################ + +#= +function jacobian( + bz::Meshes.BezierCurve, + ts::V, + diff_method::DifferentiationMethod = FiniteDifference() +) where {V <: Union{AbstractVector, Tuple}} + t = only(ts) + # Parameter t restricted to domain [0,1] by definition + if t < 0 || t > 1 + throw(DomainError(t, "b(t) is not defined for t outside [0, 1].")) + end + + # Aliases + P = bz.controls + N = Meshes.degree(bz) + + # Ensure that this implementation is tractible: limited by ability to calculate + # binomial(N, N/2) without overflow. It's possible to extend this range by + # converting N to a BigInt, but this results in always returning BigFloat's. + N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.") + + # Generator for Bernstein polynomial functions + B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i) + + # Derivative = N Σ_{i=0}^{N-1} sigma(i) + # P indices adjusted for Julia 1-based array indexing + sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1]) + derivative = N .* sum(sigma, 0:(N - 1)) + + return (derivative,) +end +=# \ No newline at end of file From e97d294dd32cb91401c2477d06aca98298ba8de4 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 27 Oct 2024 13:43:57 -0400 Subject: [PATCH 28/89] Bump CI --- src/specializations/BezierCurve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index d393797a..4f6c4f76 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -115,4 +115,4 @@ function jacobian( return (derivative,) end -=# \ No newline at end of file +=# From 7afb39c472d1c6fe76810eafb17bd0225a55e726 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 27 Oct 2024 14:29:04 -0400 Subject: [PATCH 29/89] =?UTF-8?q?Add=20B=C3=A9zierCurve=20to=20benchmarks?= =?UTF-8?q?=20(#124)=20(#126)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add Bézier entry * Missing comma, added newline * Make item implicitly Unitful * Make test item explicit to avoid reference errors * Formatting * Reformat range * Apply format suggestion --------- Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- benchmark/benchmarks.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index c438e3e2..2ec3cf17 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -19,6 +19,8 @@ rules = ( (name = "HAdaptiveCubature", rule = HAdaptiveCubature(), N = 500) ) geometries = ( + (name = "Meshes.BezierCurve", + item = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)])), (name = "Meshes.Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))), (name = "Meshes.Sphere", item = Sphere(Point(0, 0, 0), 1.0)) ) @@ -35,7 +37,7 @@ end # Differentials ############################################################################################ -sphere = geometries[2].item +sphere = Sphere(Point(0, 0, 0), 1.0) differential = MeshIntegrals.differential SUITE["Differentials"] = let s = BenchmarkGroup() From eba102b4ff4c33f968459c4416e94c49bc5d558e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 28 Oct 2024 00:18:55 -0400 Subject: [PATCH 30/89] Reduce scald --- benchmark/benchmarks.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 2ec3cf17..d830f578 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -28,7 +28,7 @@ geometries = ( SUITE["Integrals"] = let s = BenchmarkGroup() for (int, rule, geometry) in Iterators.product(integrands, rules, geometries) n1, n2, N = geometry.name, "$(int.name) $(rule.name)", rule.N - s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule) evals=N + s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule) end s end @@ -41,8 +41,8 @@ sphere = Sphere(Point(0, 0, 0), 1.0) differential = MeshIntegrals.differential SUITE["Differentials"] = let s = BenchmarkGroup() - s["Jacobian"] = @benchmarkable jacobian($sphere, $(0.5, 0.5)) evals=1000 - s["Differential"] = @benchmarkable differential($sphere, $(0.5, 0.5)) evals=1000 + s["Jacobian"] = @benchmarkable jacobian($sphere, $(0.5, 0.5)) evals=10 + s["Differential"] = @benchmarkable differential($sphere, $(0.5, 0.5)) evals=10 s end From f5a75b8c5db7576334604fc068e5c8ec23a55a25 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 28 Oct 2024 00:29:44 -0400 Subject: [PATCH 31/89] Activate specialized Jacobian --- src/specializations/BezierCurve.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 4f6c4f76..36436a5f 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -84,7 +84,6 @@ end # jacobian ################################################################################ -#= function jacobian( bz::Meshes.BezierCurve, ts::V, @@ -115,4 +114,3 @@ function jacobian( return (derivative,) end -=# From 3c93aa557b4f5635dd1dbf712b0b38db4c448f50 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 28 Oct 2024 21:13:50 -0400 Subject: [PATCH 32/89] Define an Analytic method and update jacobian of Bezier to it --- src/differentiation.jl | 2 ++ src/specializations/BezierCurve.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index c8554188..a6183c76 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -17,6 +17,8 @@ end # If ε not specified, default to 1e-6 FiniteDifference() = FiniteDifference(1e-6) +struct Analytic <: DifferentiationMethod end + # struct AutoEnzyme <: DifferentiationMethod end # struct AutoZygote <: DifferentiationMethod end diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 36436a5f..60084fe7 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -87,7 +87,7 @@ end function jacobian( bz::Meshes.BezierCurve, ts::V, - diff_method::DifferentiationMethod = FiniteDifference() + diff_method::Analytic ) where {V <: Union{AbstractVector, Tuple}} t = only(ts) # Parameter t restricted to domain [0,1] by definition From 823b56b92d97ca6711e7a674e6f0a1c087594bcd Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 28 Oct 2024 22:29:45 -0400 Subject: [PATCH 33/89] Typo [skip CI] --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index a6183c76..45326788 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -17,7 +17,7 @@ end # If ε not specified, default to 1e-6 FiniteDifference() = FiniteDifference(1e-6) -struct Analytic <: DifferentiationMethod end +struct Analytical <: DifferentiationMethod end # struct AutoEnzyme <: DifferentiationMethod end From f5c12450ed305d7af5601219771c4986d242c41e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 28 Oct 2024 22:30:09 -0400 Subject: [PATCH 34/89] Typo --- src/specializations/BezierCurve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 60084fe7..372ab742 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -87,7 +87,7 @@ end function jacobian( bz::Meshes.BezierCurve, ts::V, - diff_method::Analytic + diff_method::Analytical ) where {V <: Union{AbstractVector, Tuple}} t = only(ts) # Parameter t restricted to domain [0,1] by definition From d21b65acf8245d2f964c2e09b48128929c5345b1 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 28 Oct 2024 22:41:05 -0400 Subject: [PATCH 35/89] Update src/specializations/BezierCurve.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/specializations/BezierCurve.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 372ab742..e9ba9d6c 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -85,9 +85,9 @@ end ################################################################################ function jacobian( - bz::Meshes.BezierCurve, - ts::V, - diff_method::Analytical + bz::Meshes.BezierCurve, + ts::V, + diff_method::Analytical ) where {V <: Union{AbstractVector, Tuple}} t = only(ts) # Parameter t restricted to domain [0,1] by definition From 19e04661efa9969d496f192ef847a231d3bc3475 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 07:11:28 -0400 Subject: [PATCH 36/89] Implement has_analytic --- src/differentiation.jl | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 45326788..8aa7c910 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -1,5 +1,5 @@ ################################################################################ -# JacobianMethods +# DifferentiationMethods ################################################################################ abstract type DifferentiationMethod end @@ -19,6 +19,8 @@ FiniteDifference() = FiniteDifference(1e-6) struct Analytical <: DifferentiationMethod end +has_analytic(Type{G}) where {G <: Geometry} = false + # struct AutoEnzyme <: DifferentiationMethod end # struct AutoZygote <: DifferentiationMethod end @@ -28,23 +30,26 @@ struct Analytical <: DifferentiationMethod end ################################################################################ """ - jacobian(geometry, ts, diff_method=FiniteDifference()) + jacobian(geometry, ts[, diff_method]) -Calculate the Jacobian of a geometry's parametric function at some point `ts` -using a particular differentiation method `diff_method`. +Calculate the Jacobian of a `geometry`'s parametric function at some point `ts`. +Optionally, direct the use of a particular `differentiation method`; by default +use analytic solutions where possible and finite difference approximations +otherwise. # Arguments - `geometry`: some `Meshes.Geometry` of N parametric dimensions - `ts`: a parametric point specified as a vector or tuple of length N -- `diff_method`: the desired `DifferentiationMethod` +- `diff_method`: the desired `DifferentiationMethod` to use """ function jacobian end function jacobian( - geometry::Geometry, + geometry::G, ts::V -) where {V <: Union{AbstractVector, Tuple}} - return jacobian(geometry, ts, FiniteDifference()) +) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} + diff_method = has_analytic(G) ? Analytic() : FiniteDifference() + return jacobian(geometry, ts, diff_method) end function jacobian( @@ -87,10 +92,17 @@ end ################################################################################ """ - differential(geometry, ts, diff_method=FiniteDifference()) + differential(geometry, ts[, diff_method]) Calculate the differential element (length, area, volume, etc) of the parametric -function for `geometry` at arguments `ts`. +function for `geometry` at arguments `ts`. Optionally, direct the use of a +particular `differentiation method`; by default use analytic solutions where +possible and finite difference approximations otherwise. + +# Arguments +- `geometry`: some `Meshes.Geometry` of N parametric dimensions +- `ts`: a parametric point specified as a vector or tuple of length N +- `diff_method`: the desired `DifferentiationMethod` to use """ function differential( geometry::Geometry, From 60db456835960fc4b3ff70e096d88e4bef3e4a2f Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 07:15:57 -0400 Subject: [PATCH 37/89] Bugfix --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 8aa7c910..091a1f76 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -19,7 +19,7 @@ FiniteDifference() = FiniteDifference(1e-6) struct Analytical <: DifferentiationMethod end -has_analytic(Type{G}) where {G <: Geometry} = false +has_analytic(::Type{G}) where {G <: Geometry} = false # struct AutoEnzyme <: DifferentiationMethod end From 136ff9d9683af3acd1bd1d313d3e15e47bb5af29 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 07:17:41 -0400 Subject: [PATCH 38/89] Enable analytical for BezierCurve [skip CI] --- src/specializations/BezierCurve.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index e9ba9d6c..8b6ca06a 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -114,3 +114,5 @@ function jacobian( return (derivative,) end + +has_analytical(::BezierCurve) = true From 88b1149ea148fcb3b87344a22a596f1b40acabbf Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 07:24:39 -0400 Subject: [PATCH 39/89] Implement default_method --- src/differentiation.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 091a1f76..09647037 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -19,7 +19,8 @@ FiniteDifference() = FiniteDifference(1e-6) struct Analytical <: DifferentiationMethod end -has_analytic(::Type{G}) where {G <: Geometry} = false +has_analytical(::Type{G}) where {G <: Geometry} = false +default_method(::Type{G}) where {G <: Geometry} = has_analytical(G) ? Analytic() : FiniteDifference() # struct AutoEnzyme <: DifferentiationMethod end @@ -48,8 +49,7 @@ function jacobian( geometry::G, ts::V ) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} - diff_method = has_analytic(G) ? Analytic() : FiniteDifference() - return jacobian(geometry, ts, diff_method) + return jacobian(geometry, ts, default_method(G)) end function jacobian( @@ -107,8 +107,8 @@ possible and finite difference approximations otherwise. function differential( geometry::Geometry, ts::V, - diff_method::DifferentiationMethod = FiniteDifference() -) where {V <: Union{AbstractVector, Tuple}} + diff_method::DifferentiationMethod = default_method(G) +) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} # Calculate the Jacobian, convert Vec -> KVector J = jacobian(geometry, ts, diff_method) J_kvecs = Iterators.map(_kvector, J) From 0da8b4d4b582e26767c28e40fc2616ecbf7cf81f Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 07:42:27 -0400 Subject: [PATCH 40/89] Fix arg definition --- src/specializations/BezierCurve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 8b6ca06a..e428b9dd 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -115,4 +115,4 @@ function jacobian( return (derivative,) end -has_analytical(::BezierCurve) = true +has_analytical(::Type{Meshes.BezierCurve}) = true From 74cb3e1e060ead4a6c1c09a3a934c08096362e5f Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 07:43:27 -0400 Subject: [PATCH 41/89] Update arg type --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 09647037..7a0b217e 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -105,7 +105,7 @@ possible and finite difference approximations otherwise. - `diff_method`: the desired `DifferentiationMethod` to use """ function differential( - geometry::Geometry, + geometry::G, ts::V, diff_method::DifferentiationMethod = default_method(G) ) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} From 617a1b1a09f0adf3c2fe69e1a4edb0d3026f194d Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 19:52:14 -0400 Subject: [PATCH 42/89] Add methods on instances and types --- src/differentiation.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 7a0b217e..341bcaeb 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -19,12 +19,22 @@ FiniteDifference() = FiniteDifference(1e-6) struct Analytical <: DifferentiationMethod end +# Return whether a geometry type has Analytical() methods defined has_analytical(::Type{G}) where {G <: Geometry} = false -default_method(::Type{G}) where {G <: Geometry} = has_analytical(G) ? Analytic() : FiniteDifference() +has_analytical(g::G) where {G <: Geometry} = has_analytical(G) -# struct AutoEnzyme <: DifferentiationMethod end +# Return the default DifferentiationMethod instance for a particular geometry +function default_method( + g::Type{G} +) where {G <: Geometry} + return has_analytical(G) ? Analytic() : FiniteDifference() +end + +default_method(g::G) where {G <: Geometry} = default_method(G) -# struct AutoZygote <: DifferentiationMethod end +# Future Support: (maybe define in package extensions?) +# struct AutoEnzyme <: DifferentiationMethod end +# struct AutoZygote <: DifferentiationMethod end ################################################################################ # Jacobian From 4117374dc305a0a85fe8ed505e07ea3ff1e89ebe Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 20:05:38 -0400 Subject: [PATCH 43/89] Add tests for DifferentiationMethods --- test/utils.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/utils.jl b/test/utils.jl index 45453a50..7478ad36 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -8,4 +8,16 @@ # _units p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") @test MeshIntegrals._units(p) == u"m" + + @testset "DifferentiationMethod's" begin + bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) + sphere = Sphere(Point(0, 0, 0), 1.0) + + @test has_analytical(Meshes.BezierCurve) == true + @test has_analytical(bezier) == true + @test has_analytical(Meshes.Sphere) == false + @test has_analytical(bezier) == false + + @test FiniteDifference().ε ≈ 1e-6 + end end From bda9f3a07faa5ae83b56a89047c04dad7e3a2f3d Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 20:52:23 -0400 Subject: [PATCH 44/89] Use default_method in _integral methods --- src/integral.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index fa5032f9..0e0dba84 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - integral(f, geometry[, rule]; diff_method=FiniteDifference(), FP=Float64) + integral(f, geometry[, rule]; diff_method=default_method(geometry), FP=Float64) Numerically integrate a given function `f(::Point)` over the domain defined by a `geometry` using a particular numerical integration `rule` with floating point @@ -62,7 +62,7 @@ function _integral( f, geometry, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(geometry), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) @@ -88,7 +88,7 @@ function _integral( f, geometry, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(geometry), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) @@ -116,7 +116,7 @@ function _integral_gk_1d( f, geometry, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(geometry), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method) @@ -127,7 +127,7 @@ function _integral_gk_2d( f, geometry2d, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(geometry2d), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method) From 30d3f0b65233cd7aad938e76cf2e91dd8edaa165 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 20:52:45 -0400 Subject: [PATCH 45/89] Use default_method in Bezier methods, improve docstring --- src/specializations/BezierCurve.jl | 31 +++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index e428b9dd..8259ff3c 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -14,20 +14,29 @@ ################################################################################ """ integral(f, curve::BezierCurve, rule = GaussKronrod(); - diff_method=FiniteDifference(), FP=Float64, 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. + diff_method=Analytical(), FP=Float64, alg=Meshes.Horner()) + +Like [`integral`](@ref) but integrates along the domain defined a `curve`. + +# Arguments +- `f`: an integrand function with a method `f(::Meshes.Point)` +- `curve`: a `BezierCurve` that defines the integration domain +- `rule`: optionally, the `IntegrationRule` used for integration (by default +`GaussKronrod()` in 1D and `HAdaptiveCubature()` else) + +# Keyword Arguments +- `diff_method::DifferentiationMethod = Analytical()`: the method to use for +calculating Jacobians that are used to calculate differential elements +- `FP = Float64`: the floating point precision desired +- `alg = Meshes.Horner()`: the method to use for parameterizing `curve`. Alternatively, +`alg=Meshes.DeCasteljau()` can be specified for increased accuracy, but comes with a +steep performance cost, especially for curves with a large number of control points. """ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(curve), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} @@ -47,7 +56,7 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(curve), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} @@ -60,7 +69,7 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::DM = default_method(curve), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} From a068fd414dd8dc9922c6277939200c07939333da Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 20:55:19 -0400 Subject: [PATCH 46/89] Improve docstrings --- src/integral.jl | 12 +++++------- src/specializations/BezierCurve.jl | 3 +-- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 0e0dba84..e0ff74fd 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - integral(f, geometry[, rule]; diff_method=default_method(geometry), FP=Float64) + integral(f, geometry[, rule]; diff_method, FP) Numerically integrate a given function `f(::Point)` over the domain defined by a `geometry` using a particular numerical integration `rule` with floating point @@ -14,13 +14,11 @@ precision of type `FP`. - `geometry`: some `Meshes.Geometry` that defines the integration domain - `rule`: optionally, the `IntegrationRule` used for integration (by default `GaussKronrod()` in 1D and `HAdaptiveCubature()` else) -- `FP`: optionally, the floating point precision desired (`Float64` by default) -- `diff_method`: optionally, use a particular `DifferentiationMethod` for -calculating differential elements -Note that reducing `FP` below `Float64` will incur some loss of precision. By -contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision -(at the expense of longer runtimes). +# Keyword Arguments +- `diff_method::DifferentiationMethod = default_method(geometry)`: the method to +use for calculating Jacobians that are used to calculate differential elements +- `FP = Float64`: the floating point precision desired. """ function integral end diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 8259ff3c..b1a9dec5 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -21,8 +21,7 @@ Like [`integral`](@ref) but integrates along the domain defined a `curve`. # Arguments - `f`: an integrand function with a method `f(::Meshes.Point)` - `curve`: a `BezierCurve` that defines the integration domain -- `rule`: optionally, the `IntegrationRule` used for integration (by default -`GaussKronrod()` in 1D and `HAdaptiveCubature()` else) +- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration # Keyword Arguments - `diff_method::DifferentiationMethod = Analytical()`: the method to use for From 5b8b28313d4b9645e1b450824e387dbbfb57239e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 21:49:05 -0400 Subject: [PATCH 47/89] Explicit namespace for test of non-exported function --- test/utils.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index 7478ad36..c2365852 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -9,14 +9,14 @@ p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") @test MeshIntegrals._units(p) == u"m" - @testset "DifferentiationMethod's" begin + @testset "DifferentiationMethod" begin bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) sphere = Sphere(Point(0, 0, 0), 1.0) - @test has_analytical(Meshes.BezierCurve) == true - @test has_analytical(bezier) == true - @test has_analytical(Meshes.Sphere) == false - @test has_analytical(bezier) == false + @test Meshes.has_analytical(Meshes.BezierCurve) == true + @test Meshes.has_analytical(bezier) == true + @test Meshes.has_analytical(Meshes.Sphere) == false + @test Meshes.has_analytical(bezier) == false @test FiniteDifference().ε ≈ 1e-6 end From 42620d4e432a152314f7d6bd3a7f991ab9560016 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 22:20:06 -0400 Subject: [PATCH 48/89] Fix typo --- test/utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index c2365852..1f7f3b8f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -13,10 +13,10 @@ bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) sphere = Sphere(Point(0, 0, 0), 1.0) - @test Meshes.has_analytical(Meshes.BezierCurve) == true - @test Meshes.has_analytical(bezier) == true - @test Meshes.has_analytical(Meshes.Sphere) == false - @test Meshes.has_analytical(bezier) == false + @test MeshIntegrals.has_analytical(Meshes.BezierCurve) == true + @test MeshIntegrals.has_analytical(bezier) == true + @test MeshIntegrals.has_analytical(Meshes.Sphere) == false + @test MeshIntegrals.has_analytical(bezier) == false @test FiniteDifference().ε ≈ 1e-6 end From 969ac40e8587a2a38c93418e86e92eb62359cb93 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 22:35:50 -0400 Subject: [PATCH 49/89] Fix copy/paste error --- test/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utils.jl b/test/utils.jl index 1f7f3b8f..68461552 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -16,7 +16,7 @@ @test MeshIntegrals.has_analytical(Meshes.BezierCurve) == true @test MeshIntegrals.has_analytical(bezier) == true @test MeshIntegrals.has_analytical(Meshes.Sphere) == false - @test MeshIntegrals.has_analytical(bezier) == false + @test MeshIntegrals.has_analytical(sphere) == false @test FiniteDifference().ε ≈ 1e-6 end From f22e165378d7cf142c8c4ea478a3ec5389dc51e0 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 22:39:45 -0400 Subject: [PATCH 50/89] Add tests for default_method --- test/utils.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/utils.jl b/test/utils.jl index 68461552..cbc91db5 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -18,6 +18,11 @@ @test MeshIntegrals.has_analytical(Meshes.Sphere) == false @test MeshIntegrals.has_analytical(sphere) == false + @test MeshIntegrals.default_method(Meshes.BezierCurve) isa Analytic() + @test MeshIntegrals.default_method(bezier) isa Analytic() + @test MeshIntegrals.default_method(Meshes.Sphere) isa FiniteDifference() + @test MeshIntegrals.default_method(sphere) isa FiniteDifference() + @test FiniteDifference().ε ≈ 1e-6 end end From 2699322bc39c515fe6e30f8ad485886eb2c07f90 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 23:04:34 -0400 Subject: [PATCH 51/89] Type vs instance --- test/utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index cbc91db5..a644286f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -18,10 +18,10 @@ @test MeshIntegrals.has_analytical(Meshes.Sphere) == false @test MeshIntegrals.has_analytical(sphere) == false - @test MeshIntegrals.default_method(Meshes.BezierCurve) isa Analytic() - @test MeshIntegrals.default_method(bezier) isa Analytic() - @test MeshIntegrals.default_method(Meshes.Sphere) isa FiniteDifference() - @test MeshIntegrals.default_method(sphere) isa FiniteDifference() + @test MeshIntegrals.default_method(Meshes.BezierCurve) isa Analytic + @test MeshIntegrals.default_method(bezier) isa Analytic + @test MeshIntegrals.default_method(Meshes.Sphere) isa FiniteDifference + @test MeshIntegrals.default_method(sphere) isa FiniteDifference @test FiniteDifference().ε ≈ 1e-6 end From e8ab198ba1f6ef9145ce51ee260455f33da0cc78 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 23:31:13 -0400 Subject: [PATCH 52/89] Export Analytical --- src/MeshIntegrals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index cbd6c158..03e4dfdb 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -12,7 +12,7 @@ import Unitful include("utils.jl") include("differentiation.jl") -export DifferentiationMethod, FiniteDifference, jacobian +export DifferentiationMethod, Analytical, FiniteDifference, jacobian include("integration_rules.jl") export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature From 695d42526a4d65f24d621198ebe9f9222292d47b Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 23:31:41 -0400 Subject: [PATCH 53/89] Typo --- test/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index a644286f..ece30cc0 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -18,8 +18,8 @@ @test MeshIntegrals.has_analytical(Meshes.Sphere) == false @test MeshIntegrals.has_analytical(sphere) == false - @test MeshIntegrals.default_method(Meshes.BezierCurve) isa Analytic - @test MeshIntegrals.default_method(bezier) isa Analytic + @test MeshIntegrals.default_method(Meshes.BezierCurve) isa Analytical + @test MeshIntegrals.default_method(bezier) isa Analytical @test MeshIntegrals.default_method(Meshes.Sphere) isa FiniteDifference @test MeshIntegrals.default_method(sphere) isa FiniteDifference From b8e657e3792e2a6615b5343b90bd1686513e7679 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 23:43:55 -0400 Subject: [PATCH 54/89] Typo --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 341bcaeb..cdf483c1 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -27,7 +27,7 @@ has_analytical(g::G) where {G <: Geometry} = has_analytical(G) function default_method( g::Type{G} ) where {G <: Geometry} - return has_analytical(G) ? Analytic() : FiniteDifference() + return has_analytical(G) ? Analytical() : FiniteDifference() end default_method(g::G) where {G <: Geometry} = default_method(G) From 655555ae3a12326246ea2fc104d689d3a648c4b8 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Tue, 29 Oct 2024 23:52:32 -0400 Subject: [PATCH 55/89] Parameterize Bezier type --- src/specializations/BezierCurve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index b1a9dec5..b6cef31f 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -123,4 +123,4 @@ function jacobian( return (derivative,) end -has_analytical(::Type{Meshes.BezierCurve}) = true +has_analytical(::Type{BC}) where {BC <: Meshes.BezierCurve} = true From 8da09fda4d13f1df03eb7e1cb33d4a5b4583d496 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 30 Oct 2024 00:22:23 -0400 Subject: [PATCH 56/89] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/differentiation.jl | 2 +- test/utils.jl | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index cdf483c1..12dde5ba 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -25,7 +25,7 @@ has_analytical(g::G) where {G <: Geometry} = has_analytical(G) # Return the default DifferentiationMethod instance for a particular geometry function default_method( - g::Type{G} + g::Type{G} ) where {G <: Geometry} return has_analytical(G) ? Analytical() : FiniteDifference() end diff --git a/test/utils.jl b/test/utils.jl index ece30cc0..17d8ea0e 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -12,7 +12,6 @@ @testset "DifferentiationMethod" begin bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) sphere = Sphere(Point(0, 0, 0), 1.0) - @test MeshIntegrals.has_analytical(Meshes.BezierCurve) == true @test MeshIntegrals.has_analytical(bezier) == true @test MeshIntegrals.has_analytical(Meshes.Sphere) == false From 8f528821eb708114b8acb1f365860aafb0795c1c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 30 Oct 2024 14:21:32 -0400 Subject: [PATCH 57/89] Redefine Line, Plane, and Ray specializations as Analytical --- src/specializations/Line.jl | 8 +++++--- src/specializations/Plane.jl | 8 +++++--- src/specializations/Ray.jl | 8 +++++--- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 25617474..480ff484 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -11,7 +11,7 @@ function integral( f::F, line::Meshes.Line, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] @@ -34,7 +34,7 @@ function integral( f::F, line::Meshes.Line, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Line s.t. line(t) is distance t from origin point @@ -50,7 +50,7 @@ function integral( f::F, line::Meshes.Line, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Line s.t. line(t) is distance t from origin point @@ -76,3 +76,5 @@ function integral( # Reapply units return value .* integrandunits end + +has_analytical(::Type{T}) where {T <: Meshes.Line} = true diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 79039868..2323576a 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -12,7 +12,7 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² @@ -41,7 +41,7 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Plane's orthogonal vectors @@ -60,7 +60,7 @@ function integral( f::F, plane::Meshes.Plane, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Plane's orthogonal vectors @@ -90,3 +90,5 @@ function integral( # Reapply units return value .* integrandunits end + +has_analytical(::Type{T}) where {T <: Meshes.Plane} = true diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index e17f64a6..7c9ae7c0 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -12,7 +12,7 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] @@ -39,7 +39,7 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Ray s.t. ray(t) is distance t from origin point @@ -55,7 +55,7 @@ function integral( f::F, ray::Meshes.Ray, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Normalize the Ray s.t. ray(t) is distance t from origin point @@ -81,3 +81,5 @@ function integral( # Reapply units return value .* integrandunits end + +has_analytical(::Type{T}) where {T <: Meshes.Ray} = true From 678e8cc2d321fb4942ca3e9fff828385eb64fb87 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 30 Oct 2024 14:23:03 -0400 Subject: [PATCH 58/89] Redefine Triangle and Tetrahedron specializations as Analytical --- src/specializations/Tetrahedron.jl | 8 +++++--- src/specializations/Triangle.jl | 8 +++++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index 4c9c52d3..58523c9f 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -13,7 +13,7 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _error_unsupported_combination("Tetrahedron", "GaussLegendre") @@ -23,7 +23,7 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} nil = zero(FP) @@ -40,8 +40,10 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature") end + +has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index acca16b2..1f74cea8 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -21,7 +21,7 @@ function integral( f::F, triangle::Meshes.Ngon{3}, rule::GaussLegendre; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 @@ -64,7 +64,7 @@ function integral( f::F, triangle::Meshes.Ngon{3}, rule::GaussKronrod; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0) @@ -87,7 +87,7 @@ function integral( f::F, triangle::Meshes.Ngon{3}, rule::HAdaptiveCubature; - diff_method::DM = FiniteDifference(), + diff_method::Analytical, FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} # Integrate the Barycentric triangle by transforming it into polar coordinates @@ -107,3 +107,5 @@ function integral( # Apply a linear domain-correction factor 0.5 ↦ area(triangle) return 2 * Meshes.area(triangle) .* ∫ end + +has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true From 1560ebfc68042f4f1bafcd3a320bc9d288e9597a Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 08:33:33 -0400 Subject: [PATCH 59/89] Try new solution (cant dispatch on kwarg types) --- src/specializations/Line.jl | 12 +++++++++--- src/utils.jl | 6 ++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 480ff484..8659ffd5 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -11,9 +11,11 @@ function integral( f::F, line::Meshes.Line, rule::GaussLegendre; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Line, diff_method) + # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] xs, ws = _gausslegendre(FP, rule.n) @@ -34,9 +36,11 @@ function integral( f::F, line::Meshes.Line, rule::GaussKronrod; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Line, diff_method) + # Normalize the Line s.t. line(t) is distance t from origin point line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) @@ -50,9 +54,11 @@ function integral( f::F, line::Meshes.Line, rule::HAdaptiveCubature; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Line, diff_method) + # Normalize the Line s.t. line(t) is distance t from origin point line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) diff --git a/src/utils.jl b/src/utils.jl index 14de6e6b..5243eec5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,6 +14,12 @@ function _error_unsupported_combination(geometry, rule) throw(ArgumentError(msg)) end +function _guarantee_analytical(G, diff_method) + if diff_method != Analytical() + ArgumentError("Geometry type $G requires kwarg diff_method = Analytical()") + end +end + ################################################################################ # CliffordNumbers and Units ################################################################################ From fee4d3503f739157bb5a35488aeb1868482a81e7 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 12:55:53 -0400 Subject: [PATCH 60/89] Add jacobian code sections --- src/specializations/Line.jl | 4 ++++ src/specializations/Plane.jl | 4 ++++ src/specializations/Ray.jl | 4 ++++ src/specializations/Tetrahedron.jl | 4 ++++ src/specializations/Triangle.jl | 4 ++++ 5 files changed, 20 insertions(+) diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 8659ffd5..402b812e 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -83,4 +83,8 @@ function integral( return value .* integrandunits end +################################################################################ +# jacobian +################################################################################ + has_analytical(::Type{T}) where {T <: Meshes.Line} = true diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 2323576a..7bc773b0 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -91,4 +91,8 @@ function integral( return value .* integrandunits end +################################################################################ +# jacobian +################################################################################ + has_analytical(::Type{T}) where {T <: Meshes.Plane} = true diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index 7c9ae7c0..3ae4def3 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -82,4 +82,8 @@ function integral( return value .* integrandunits end +################################################################################ +# jacobian +################################################################################ + has_analytical(::Type{T}) where {T <: Meshes.Ray} = true diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index 58523c9f..d86c1a64 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -46,4 +46,8 @@ function integral( _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature") end +################################################################################ +# jacobian +################################################################################ + has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index 1f74cea8..5a57f334 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -108,4 +108,8 @@ function integral( return 2 * Meshes.area(triangle) .* ∫ end +################################################################################ +# jacobian +################################################################################ + has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true From 73f0b9abecd48e2cfccbeb52a3f0b54c2c9ffa84 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 12:59:38 -0400 Subject: [PATCH 61/89] Add guarantees --- src/specializations/BezierCurve.jl | 2 +- src/specializations/Plane.jl | 12 +++++++++--- src/specializations/Ray.jl | 12 +++++++++--- src/specializations/Tetrahedron.jl | 8 +++++--- src/specializations/Triangle.jl | 18 ++++++++++++------ 5 files changed, 36 insertions(+), 16 deletions(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index b6cef31f..ac411b89 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -123,4 +123,4 @@ function jacobian( return (derivative,) end -has_analytical(::Type{BC}) where {BC <: Meshes.BezierCurve} = true +has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 7bc773b0..b46a4b57 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -12,9 +12,11 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussLegendre; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Plane, diff_method) + # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² xs, ws = _gausslegendre(FP, rule.n) wws = Iterators.product(ws, ws) @@ -41,9 +43,11 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussKronrod; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Plane, diff_method) + # Normalize the Plane's orthogonal vectors uu = Meshes.unormalize(plane.u) uv = Meshes.unormalize(plane.v) @@ -60,9 +64,11 @@ function integral( f::F, plane::Meshes.Plane, rule::HAdaptiveCubature; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Plane, diff_method) + # Normalize the Plane's orthogonal vectors uu = Meshes.unormalize(plane.u) uv = Meshes.unormalize(plane.v) diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index 3ae4def3..3ded2a90 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -12,9 +12,11 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussLegendre; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Ray, diff_method) + # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] xs, ws = _gausslegendre(FP, rule.n) @@ -39,9 +41,11 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussKronrod; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Ray, diff_method) + # Normalize the Ray s.t. ray(t) is distance t from origin point ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) @@ -55,9 +59,11 @@ function integral( f::F, ray::Meshes.Ray, rule::HAdaptiveCubature; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Ray, diff_method) + # Normalize the Ray s.t. ray(t) is distance t from origin point ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index d86c1a64..50457525 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -13,7 +13,7 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussLegendre; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _error_unsupported_combination("Tetrahedron", "GaussLegendre") @@ -23,9 +23,11 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussKronrod; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Tetrahedron, diff_method) + nil = zero(FP) ∫uvw(u, v, w) = f(tetrahedron(u, v, w)) ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), nil, FP(1 - v - w); rule.kwargs...)[1] @@ -40,7 +42,7 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::HAdaptiveCubature; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature") diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index 5a57f334..2199d853 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -19,11 +19,13 @@ triangle. """ function integral( f::F, - triangle::Meshes.Ngon{3}, + triangle::Meshes.Triangle, rule::GaussLegendre; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Triangle, diff_method) + # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 xs, ws = _gausslegendre(FP, rule.n) wws = Iterators.product(ws, ws) @@ -62,11 +64,13 @@ Gauss-Kronrod quadrature rules along each barycentric dimension of the triangle. """ function integral( f::F, - triangle::Meshes.Ngon{3}, + triangle::Meshes.Triangle, rule::GaussKronrod; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Triangle, diff_method) + # 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 ∫u(u) = QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1] @@ -85,11 +89,13 @@ an h-adaptive cubature rule. """ function integral( f::F, - triangle::Meshes.Ngon{3}, + triangle::Meshes.Triangle, rule::HAdaptiveCubature; - diff_method::Analytical, + diff_method::DM = Analytical(), FP::Type{T} = Float64 ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} + _guarantee_analytical(Meshes.Triangle, diff_method) + # Integrate the Barycentric triangle by transforming it into polar coordinates # with a modified radius # R = r ( sinφ + cosφ ) From 77c36978fa13fb9ead1e5b2919e67b7fc1365631 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 19:40:36 -0400 Subject: [PATCH 62/89] Move new utilities to utils and rename with prefix underscore --- src/differentiation.jl | 13 ------------- src/integral.jl | 10 +++++----- src/specializations/BezierCurve.jl | 8 ++++---- src/specializations/Line.jl | 2 +- src/specializations/Plane.jl | 2 +- src/specializations/Ray.jl | 2 +- src/specializations/Tetrahedron.jl | 2 +- src/specializations/Triangle.jl | 2 +- src/utils.jl | 19 +++++++++++++++++++ 9 files changed, 33 insertions(+), 27 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 12dde5ba..0e377cd3 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -19,19 +19,6 @@ FiniteDifference() = FiniteDifference(1e-6) struct Analytical <: DifferentiationMethod end -# Return whether a geometry type has Analytical() methods defined -has_analytical(::Type{G}) where {G <: Geometry} = false -has_analytical(g::G) where {G <: Geometry} = has_analytical(G) - -# Return the default DifferentiationMethod instance for a particular geometry -function default_method( - g::Type{G} -) where {G <: Geometry} - return has_analytical(G) ? Analytical() : FiniteDifference() -end - -default_method(g::G) where {G <: Geometry} = default_method(G) - # Future Support: (maybe define in package extensions?) # struct AutoEnzyme <: DifferentiationMethod end # struct AutoZygote <: DifferentiationMethod end diff --git a/src/integral.jl b/src/integral.jl index e0ff74fd..30f3bc85 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -16,7 +16,7 @@ precision of type `FP`. `GaussKronrod()` in 1D and `HAdaptiveCubature()` else) # Keyword Arguments -- `diff_method::DifferentiationMethod = default_method(geometry)`: the method to +- `diff_method::DifferentiationMethod = _default_method(geometry)`: the method to use for calculating Jacobians that are used to calculate differential elements - `FP = Float64`: the floating point precision desired. """ @@ -60,7 +60,7 @@ function _integral( f, geometry, rule::GaussLegendre; - diff_method::DM = default_method(geometry), + diff_method::DM = _default_method(geometry), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) @@ -86,7 +86,7 @@ function _integral( f, geometry, rule::HAdaptiveCubature; - diff_method::DM = default_method(geometry), + diff_method::DM = _default_method(geometry), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) @@ -114,7 +114,7 @@ function _integral_gk_1d( f, geometry, rule::GaussKronrod; - diff_method::DM = default_method(geometry), + diff_method::DM = _default_method(geometry), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method) @@ -125,7 +125,7 @@ function _integral_gk_2d( f, geometry2d, rule::GaussKronrod; - diff_method::DM = default_method(geometry2d), + diff_method::DM = _default_method(geometry2d), FP::Type{T} = Float64 ) where {DM <: DifferentiationMethod, T <: AbstractFloat} integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index ac411b89..4de250b3 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -35,7 +35,7 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussLegendre; - diff_method::DM = default_method(curve), + diff_method::DM = _default_method(curve), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} @@ -55,7 +55,7 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussKronrod; - diff_method::DM = default_method(curve), + diff_method::DM = _default_method(curve), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} @@ -68,7 +68,7 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::HAdaptiveCubature; - diff_method::DM = default_method(curve), + diff_method::DM = _default_method(curve), FP::Type{T} = Float64, alg::Meshes.BezierEvalMethod = Meshes.Horner() ) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} @@ -123,4 +123,4 @@ function jacobian( return (derivative,) end -has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true +_has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 402b812e..655278e1 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -87,4 +87,4 @@ end # jacobian ################################################################################ -has_analytical(::Type{T}) where {T <: Meshes.Line} = true +_has_analytical(::Type{T}) where {T <: Meshes.Line} = true diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index b46a4b57..1b34405e 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -101,4 +101,4 @@ end # jacobian ################################################################################ -has_analytical(::Type{T}) where {T <: Meshes.Plane} = true +_has_analytical(::Type{T}) where {T <: Meshes.Plane} = true diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index 3ded2a90..9008b4a8 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -92,4 +92,4 @@ end # jacobian ################################################################################ -has_analytical(::Type{T}) where {T <: Meshes.Ray} = true +_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index 50457525..9fb2af9a 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -52,4 +52,4 @@ end # jacobian ################################################################################ -has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true +_has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index 2199d853..20c68ea3 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -118,4 +118,4 @@ end # jacobian ################################################################################ -has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true +_has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true diff --git a/src/utils.jl b/src/utils.jl index 5243eec5..213bbb7d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,12 +14,31 @@ function _error_unsupported_combination(geometry, rule) throw(ArgumentError(msg)) end +################################################################################ +# DifferentiationMethod +################################################################################ + +# Throw an ArgumentError if Analytical() jacobian not defined for this type function _guarantee_analytical(G, diff_method) if diff_method != Analytical() ArgumentError("Geometry type $G requires kwarg diff_method = Analytical()") end end +# Return whether a geometry type has jacobian methods defined +_has_analytical(::Type{G}) where {G <: Geometry} = false +_has_analytical(g::G) where {G <: Geometry} = has_analytical(G) + +# Return the default DifferentiationMethod instance for a particular geometry type +function _default_method( + g::Type{G} +) where {G <: Geometry} + return _has_analytical(G) ? Analytical() : FiniteDifference() +end + +# Return the default DifferentiationMethod instance for a particular geometry instance +_default_method(g::G) where {G <: Geometry} = default_method(G) + ################################################################################ # CliffordNumbers and Units ################################################################################ From 94f18f516c885391e6430890d5d673e513eb402a Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 19:52:19 -0400 Subject: [PATCH 63/89] Typo --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 213bbb7d..32d3bb72 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -37,7 +37,7 @@ function _default_method( end # Return the default DifferentiationMethod instance for a particular geometry instance -_default_method(g::G) where {G <: Geometry} = default_method(G) +_default_method(g::G) where {G <: Geometry} = _default_method(G) ################################################################################ # CliffordNumbers and Units From 1af6a60f38f79ea40c08ccb2a494b81caad22e57 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 19:52:26 -0400 Subject: [PATCH 64/89] Update tests --- test/utils.jl | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index 17d8ea0e..5706f097 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -10,17 +10,25 @@ @test MeshIntegrals._units(p) == u"m" @testset "DifferentiationMethod" begin + # Test _has_analytical of instances bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) sphere = Sphere(Point(0, 0, 0), 1.0) - @test MeshIntegrals.has_analytical(Meshes.BezierCurve) == true - @test MeshIntegrals.has_analytical(bezier) == true - @test MeshIntegrals.has_analytical(Meshes.Sphere) == false - @test MeshIntegrals.has_analytical(sphere) == false + @test MeshIntegrals._has_analytical(bezier) == true + @test MeshIntegrals._has_analytical(sphere) == false - @test MeshIntegrals.default_method(Meshes.BezierCurve) isa Analytical - @test MeshIntegrals.default_method(bezier) isa Analytical - @test MeshIntegrals.default_method(Meshes.Sphere) isa FiniteDifference - @test MeshIntegrals.default_method(sphere) isa FiniteDifference + # Test _has_analytical of types + @test MeshIntegrals._has_analytical(Meshes.BezierCurve) == true + @test MeshIntegrals._has_analytical(Meshes.Line) == true + @test MeshIntegrals._has_analytical(Meshes.Plane) == true + @test MeshIntegrals._has_analytical(Meshes.Ray) == true + @test MeshIntegrals._has_analytical(Meshes.Sphere) == false + @test MeshIntegrals._has_analytical(Meshes.Tetrahedron) == true + @test MeshIntegrals._has_analytical(Meshes.Triangle) == true + + @test MeshIntegrals._default_method(Meshes.BezierCurve) isa Analytical + @test MeshIntegrals._default_method(bezier) isa Analytical + @test MeshIntegrals._default_method(Meshes.Sphere) isa FiniteDifference + @test MeshIntegrals._default_method(sphere) isa FiniteDifference @test FiniteDifference().ε ≈ 1e-6 end From 36a0426a45e5da19cc996770ec32b9cbf4276872 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 19:53:03 -0400 Subject: [PATCH 65/89] Update name --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 32d3bb72..8404867e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -27,7 +27,7 @@ end # Return whether a geometry type has jacobian methods defined _has_analytical(::Type{G}) where {G <: Geometry} = false -_has_analytical(g::G) where {G <: Geometry} = has_analytical(G) +_has_analytical(g::G) where {G <: Geometry} = _has_analytical(G) # Return the default DifferentiationMethod instance for a particular geometry type function _default_method( From afbd927fb3b7e07709840c37f2bf038219030d74 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 20:02:12 -0400 Subject: [PATCH 66/89] Update function name --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 0e377cd3..69660163 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -46,7 +46,7 @@ function jacobian( geometry::G, ts::V ) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} - return jacobian(geometry, ts, default_method(G)) + return jacobian(geometry, ts, _default_method(G)) end function jacobian( From a36b86466d5ab51d6c52fd08eb1127d6e95fc4f3 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 20:02:47 -0400 Subject: [PATCH 67/89] Update comment --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 69660163..2d101352 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -19,7 +19,7 @@ FiniteDifference() = FiniteDifference(1e-6) struct Analytical <: DifferentiationMethod end -# Future Support: (maybe define in package extensions?) +# Future Support: # struct AutoEnzyme <: DifferentiationMethod end # struct AutoZygote <: DifferentiationMethod end From 613358958c641885cbfd44e5ea7be1ec2c8e1e70 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 20:43:01 -0400 Subject: [PATCH 68/89] Update name --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 2d101352..5c6b3cb8 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -104,7 +104,7 @@ possible and finite difference approximations otherwise. function differential( geometry::G, ts::V, - diff_method::DifferentiationMethod = default_method(G) + diff_method::DifferentiationMethod = _default_method(G) ) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} # Calculate the Jacobian, convert Vec -> KVector J = jacobian(geometry, ts, diff_method) From 06946a133a1479916b20213dabc00e73f1806a92 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 21:01:45 -0400 Subject: [PATCH 69/89] Add better docstrings --- src/differentiation.jl | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 5c6b3cb8..7cfb074c 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -2,13 +2,20 @@ # DifferentiationMethods ################################################################################ +""" + DifferentiationMethod + +A category of types used to specify the desired method for calculating derivatives. +Derivatives are used to form Jacobian matrices when calculating the differential +element size throughout the integration region. +""" abstract type DifferentiationMethod end """ FiniteDifference(ε=1e-6) -Use a finite-difference approximation method to calculate derivatives with a -step size of `ε`. +Use to specify use of a finite-difference approximation method with a step size +of `ε` for calculating derivatives. """ struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod ε::T @@ -17,6 +24,20 @@ end # If ε not specified, default to 1e-6 FiniteDifference() = FiniteDifference(1e-6) +""" + Analytical() + +Use to specify use of analytically-derived solutions for calculating derivatives. +These solutions are currently defined only for a subset of geometry types. + +# Supported Geometries: +- `BezierCurve` +- `Line` +- `Plane` +- `Ray` +- `Tetrahedron` +- `Triangle` +""" struct Analytical <: DifferentiationMethod end # Future Support: From cc31a296e0c9feab200d34fbd6530c304d703980 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 21:21:45 -0400 Subject: [PATCH 70/89] Bugfix - error not actually thrown --- src/utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 8404867e..551ac1e1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -19,12 +19,12 @@ end ################################################################################ # Throw an ArgumentError if Analytical() jacobian not defined for this type -function _guarantee_analytical(G, diff_method) - if diff_method != Analytical() - ArgumentError("Geometry type $G requires kwarg diff_method = Analytical()") - end +function _guarantee_analytical(G, diff_method::DifferentiationMethod) + throw(ArgumentError("Geometry type $G requires kwarg diff_method = Analytical()")) end +_guarantee_analytical(G, ::Analytical) = nothing + # Return whether a geometry type has jacobian methods defined _has_analytical(::Type{G}) where {G <: Geometry} = false _has_analytical(g::G) where {G <: Geometry} = _has_analytical(G) From 45b509c98e17cf916cf096f1c1fd43178825d60f Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 21:22:03 -0400 Subject: [PATCH 71/89] Split DifferentiationMethod into separate testitem --- test/utils.jl | 56 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 23 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index 5706f097..2ab4ac19 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -9,27 +9,37 @@ p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") @test MeshIntegrals._units(p) == u"m" - @testset "DifferentiationMethod" begin - # Test _has_analytical of instances - bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) - sphere = Sphere(Point(0, 0, 0), 1.0) - @test MeshIntegrals._has_analytical(bezier) == true - @test MeshIntegrals._has_analytical(sphere) == false - - # Test _has_analytical of types - @test MeshIntegrals._has_analytical(Meshes.BezierCurve) == true - @test MeshIntegrals._has_analytical(Meshes.Line) == true - @test MeshIntegrals._has_analytical(Meshes.Plane) == true - @test MeshIntegrals._has_analytical(Meshes.Ray) == true - @test MeshIntegrals._has_analytical(Meshes.Sphere) == false - @test MeshIntegrals._has_analytical(Meshes.Tetrahedron) == true - @test MeshIntegrals._has_analytical(Meshes.Triangle) == true - - @test MeshIntegrals._default_method(Meshes.BezierCurve) isa Analytical - @test MeshIntegrals._default_method(bezier) isa Analytical - @test MeshIntegrals._default_method(Meshes.Sphere) isa FiniteDifference - @test MeshIntegrals._default_method(sphere) isa FiniteDifference - - @test FiniteDifference().ε ≈ 1e-6 - end + end + +@testitem "DifferentiationMethod" setup=[Setup] begin + using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical + + # _has_analytical of instances + bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) + @test _has_analytical(bezier) == true + sphere = Sphere(Point(0, 0, 0), 1.0) + @test _has_analytical(sphere) == false + + # _has_analytical of types + @test _has_analytical(Meshes.BezierCurve) == true + @test _has_analytical(Meshes.Line) == true + @test _has_analytical(Meshes.Plane) == true + @test _has_analytical(Meshes.Ray) == true + @test _has_analytical(Meshes.Sphere) == false + @test _has_analytical(Meshes.Tetrahedron) == true + @test _has_analytical(Meshes.Triangle) == true + + # _guarantee_analytical + @test _guarantee_analytical(Meshes.Line, Analytical()) === nothing + @test_throws "Analytical" _guarantee_analytical(Meshes.Line, FiniteDifference()) + + # _default_method + @test _default_method(Meshes.BezierCurve) isa Analytical + @test _default_method(bezier) isa Analytical + @test _default_method(Meshes.Sphere) isa FiniteDifference + @test _default_method(sphere) isa FiniteDifference + + # FiniteDifference + @test FiniteDifference().ε ≈ 1e-6 +end \ No newline at end of file From e6c226f4ae05610559d7db834ebfd542544fb46e Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 21:28:09 -0400 Subject: [PATCH 72/89] Tweak include order --- src/MeshIntegrals.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 03e4dfdb..5fbac40f 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -9,11 +9,11 @@ import LinearAlgebra import QuadGK import Unitful -include("utils.jl") - include("differentiation.jl") export DifferentiationMethod, Analytical, FiniteDifference, jacobian +include("utils.jl") + include("integration_rules.jl") export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature From fb6c370e8fc6d1f567b611227a22a76df68b59d2 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 21:29:25 -0400 Subject: [PATCH 73/89] Formatting changes --- test/utils.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index 2ab4ac19..2a2d8a68 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -8,8 +8,6 @@ # _units p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") @test MeshIntegrals._units(p) == u"m" - - end @testitem "DifferentiationMethod" setup=[Setup] begin @@ -42,4 +40,4 @@ end # FiniteDifference @test FiniteDifference().ε ≈ 1e-6 -end \ No newline at end of file +end From 4db8280a470d0c3f26422b9ee02b6495065484cc Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 23:23:07 -0400 Subject: [PATCH 74/89] Remove redundant docstrings --- src/specializations/Triangle.jl | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index 20c68ea3..c63e715b 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -9,14 +9,6 @@ # documentation. ################################################################################ -""" - integral(f, triangle::Triangle, ::GaussLegendre; FP=Float64) - -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.Triangle, @@ -56,12 +48,6 @@ function integral( return FP(π / 4) * Meshes.area(triangle) .* sum(integrand, zip(wws, xxs)) end -""" - integral(f, triangle::Triangle, ::GaussKronrod; FP=Float64) - -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.Triangle, @@ -80,13 +66,6 @@ function integral( return 2 * Meshes.area(triangle) .* ∫ end -""" - integral(f, triangle::Triangle, ::GaussKronrod; FP=Float64) - -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.Triangle, From 65db4f64f366d405a0248a2e949af0886b3042a7 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Fri, 1 Nov 2024 23:28:02 -0400 Subject: [PATCH 75/89] Add docstrings for Ring and Rope, update for Bezier --- src/specializations/BezierCurve.jl | 2 +- src/specializations/Ring.jl | 18 ++++++++++++++++++ src/specializations/Rope.jl | 18 ++++++++++++++++++ 3 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 4de250b3..a6b06828 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -16,7 +16,7 @@ integral(f, curve::BezierCurve, rule = GaussKronrod(); diff_method=Analytical(), FP=Float64, alg=Meshes.Horner()) -Like [`integral`](@ref) but integrates along the domain defined a `curve`. +Like [`integral`](@ref) but integrates along the domain defined by `curve`. # Arguments - `f`: an integrand function with a method `f(::Meshes.Point)` diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl index 51cadf56..b2490ee2 100644 --- a/src/specializations/Ring.jl +++ b/src/specializations/Ring.jl @@ -8,6 +8,24 @@ # constituent Segment's, integrated separately, and then summed. ################################################################################ +""" + integral(f, ring::Ring, rule = GaussKronrod(); + diff_method=FiniteDifference(), FP=Float64) + +Like [`integral`](@ref) but integrates along the domain defined by `ring`. The +specified integration `rule` is applied independently to each segment formed by +consecutive points in the Ring. + +# Arguments +- `f`: an integrand function with a method `f(::Meshes.Point)` +- `ring`: a `Ring` that defines the integration domain +- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration + +# Keyword Arguments +- `diff_method::DifferentiationMethod = FiniteDifference()`: the method to use for +calculating Jacobians that are used to calculate differential elements +- `FP = Float64`: the floating point precision desired +""" function integral( f::F, ring::Meshes.Ring, diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl index c68e82cc..7fcefa7f 100644 --- a/src/specializations/Rope.jl +++ b/src/specializations/Rope.jl @@ -8,6 +8,24 @@ # integrated separately, and then summed. ################################################################################ +""" + integral(f, rope::Rope, rule = GaussKronrod(); + diff_method=FiniteDifference(), FP=Float64) + +Like [`integral`](@ref) but integrates along the domain defined by `rope`. The +specified integration `rule` is applied independently to each segment formed by +consecutive points in the Rope. + +# Arguments +- `f`: an integrand function with a method `f(::Meshes.Point)` +- `rope`: a `Rope` that defines the integration domain +- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration + +# Keyword Arguments +- `diff_method::DifferentiationMethod = FiniteDifference()`: the method to use for +calculating Jacobians that are used to calculate differential elements +- `FP = Float64`: the floating point precision desired +""" function integral( f::F, rope::Meshes.Rope, From 922518862cbb707a4eed000b49d9e388fa096716 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 09:15:06 -0400 Subject: [PATCH 76/89] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/differentiation.jl | 12 ++++++------ src/integral.jl | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 7cfb074c..d0941129 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -8,6 +8,8 @@ A category of types used to specify the desired method for calculating derivatives. Derivatives are used to form Jacobian matrices when calculating the differential element size throughout the integration region. + +See also [`FiniteDifference`](@ref), [`Analytical`](@ref). """ abstract type DifferentiationMethod end @@ -52,17 +54,15 @@ struct Analytical <: DifferentiationMethod end jacobian(geometry, ts[, diff_method]) Calculate the Jacobian of a `geometry`'s parametric function at some point `ts`. -Optionally, direct the use of a particular `differentiation method`; by default +Optionally, direct the use of a particular differentiation method `diff_method`; by default use analytic solutions where possible and finite difference approximations otherwise. # Arguments - `geometry`: some `Meshes.Geometry` of N parametric dimensions - `ts`: a parametric point specified as a vector or tuple of length N -- `diff_method`: the desired `DifferentiationMethod` to use +- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use """ -function jacobian end - function jacobian( geometry::G, ts::V @@ -114,13 +114,13 @@ end Calculate the differential element (length, area, volume, etc) of the parametric function for `geometry` at arguments `ts`. Optionally, direct the use of a -particular `differentiation method`; by default use analytic solutions where +particular differentiation method `diff_method`; by default use analytic solutions where possible and finite difference approximations otherwise. # Arguments - `geometry`: some `Meshes.Geometry` of N parametric dimensions - `ts`: a parametric point specified as a vector or tuple of length N -- `diff_method`: the desired `DifferentiationMethod` to use +- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use """ function differential( geometry::G, diff --git a/src/integral.jl b/src/integral.jl index 30f3bc85..f0f9389a 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - integral(f, geometry[, rule]; diff_method, FP) + integral(f, geometry[, rule]; diff_method=_default_method(geometry), FP=Float64) Numerically integrate a given function `f(::Point)` over the domain defined by a `geometry` using a particular numerical integration `rule` with floating point From 63bc35ec17d3adae50e6e452906649bed6c469b5 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 14:50:46 -0400 Subject: [PATCH 77/89] (WIP) Add a new kwarg section under Usage --- README.md | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 68cc0ce0..6bc7d144 100644 --- a/README.md +++ b/README.md @@ -20,17 +20,19 @@ These solvers have support for integrand functions that produce scalars, vectors ## Usage +### Common API + ```julia integral(f, geometry) ``` Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry`. A default integration method will be automatically selected according to the geometry: `GaussKronrod()` for 1D, and `HAdaptiveCubature()` for all others. ```julia -integral(f, geometry, rule, FP=Float64) +integral(f, geometry, rule) ``` Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. -Optionally, a fourth argument can be provided to specify the floating point precision level desired. This setting can be manipulated if your integrand function produces outputs with alternate floating point precision (e.g. `Float16`, `BigFloat`, etc) AND you'd prefer to avoid implicit type promotions. +### Aliases ```julia lineintegral(f, geometry) @@ -42,6 +44,14 @@ Alias functions are provided for convenience. These are simply wrappers for `int - `surfaceintegral` for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc) - `volumeintegral` for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc) +### Optional Keyword Arguments + +The `diff_method` keyword argument ... + +The `FP` keyword argument can be provided to specify the floating point precision level desired. This setting can be manipulated if your integrand function produces outputs with alternate floating point precision (e.g. `Float16`, `BigFloat`, etc) AND you'd prefer to avoid implicit type promotions. + + + # Demo ```julia From 0e4db4b0418a36e9c1e807330dd36b8c2fc13c45 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 21:54:10 -0400 Subject: [PATCH 78/89] Clarify Usage section --- README.md | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 6bc7d144..6c8ac1a0 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ These solvers have support for integrand functions that produce scalars, vectors ## Usage -### Common API +### Basic ```julia integral(f, geometry) @@ -32,6 +32,9 @@ integral(f, geometry, rule) ``` Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. +Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) +to provide additional control over the integration mechanics. + ### Aliases ```julia @@ -39,18 +42,10 @@ lineintegral(f, geometry) surfaceintegral(f, geometry) volumeintegral(f, geometry) ``` -Alias functions are provided for convenience. These are simply wrappers for `integral` that first validate that the provided `geometry` has the expected number of parametric/manifold dimensions. Like with `integral` in the examples above, the `rule` can also be specified as a third-argument. -- `lineintegral` for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc) -- `surfaceintegral` for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc) -- `volumeintegral` for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc) - -### Optional Keyword Arguments - -The `diff_method` keyword argument ... - -The `FP` keyword argument can be provided to specify the floating point precision level desired. This setting can be manipulated if your integrand function produces outputs with alternate floating point precision (e.g. `Float16`, `BigFloat`, etc) AND you'd prefer to avoid implicit type promotions. - - +Alias functions are provided for convenience. These are simply wrappers for `integral` that also validate that the provided `geometry` has the expected number of parametric dimensions. Like with `integral`, a `rule` can also optionally be specified as a third argument. +- `lineintegral` is used for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc) +- `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc) +- `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc) # Demo From 4930a2b75f54a0c87da0d1bd2a0808c71328172a Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 21:57:13 -0400 Subject: [PATCH 79/89] Formatting --- README.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/README.md b/README.md index 6c8ac1a0..1e51d48c 100644 --- a/README.md +++ b/README.md @@ -32,8 +32,7 @@ integral(f, geometry, rule) ``` Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. -Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) -to provide additional control over the integration mechanics. +Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics. ### Aliases From 6907489fb92ee4e4000e2b0a0515f8b688c1df97 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 22:30:38 -0400 Subject: [PATCH 80/89] Add a WIP differential forms explainer --- docs/src/differential_forms.md | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 docs/src/differential_forms.md diff --git a/docs/src/differential_forms.md b/docs/src/differential_forms.md new file mode 100644 index 00000000..b8e34ab3 --- /dev/null +++ b/docs/src/differential_forms.md @@ -0,0 +1,30 @@ +# Differential Forms (Work in Progress) + +**MeshIntegrals.jl** uses differential forms to perform numerical integration. Every supported `Meshes.Geometry` type is defined as having a parametric function that maps from a local coordinate system to every point on the geometry. For example, +```julia +a = Meshes.Point(0, 0) +b = Meshes.Point(2, 4) +segment = Meshes.Segment(a, b) +``` +defines a line segment beginning at point `a` and ending at point `b`. As a geometry with one parametric dimension (i.e. `paramdim(segment) == 1`), it can be treated as a function with one argument that returns a corresponding point on the segment. +```julia +segment(0) == a +segment(0.5) == Point(1, 2) +segment(1) == b +``` + +The general solution for integrating such a geometry in 3D space can look something like the following, where $t$ is a parametric coordinate used to generate points in the domain. + +**TODO: update all of the above for a 2D geometry (Sphere?) to make it more interesting/relevant. Something simple enough to follow but non-trivial.** + +Using differential forms, the general solution for integrating a geometry with two parametric dimensions, $t_1$ and $t_2$, is +```math +\iint f(x, y, z) ~ \text{d}t_1 \wedge \text{d}t_2 +``` + +Since `Meshes.Geometry`s parametric functions have arguments that are defined on the domain $[0,1]$, this is equivalent to +```math +\int_0^1 \int_0^1 f(x, y, z) ~ \text{d}t_1 \wedge \text{d}t_2 +``` + +For every point in the integration domain where the integrand function is evaluated, the differential element $\text{d}t_1 \wedge \text{d}t_2$ is calculated using [the Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of the parametric function. For a two-dimensional surface this Jacobian consists of two vectors, each pointing in the direction that the parametric function's output point will move by changing each input argument. The differential element, then, is simply the magnitude of the exterior product of these vectors. \ No newline at end of file From 2b462d13d4090c4970aba6a969bbc8f57d915660 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 23:34:03 -0400 Subject: [PATCH 81/89] Rename page and add some structure --- ...{differential_forms.md => how_it_works.md} | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) rename docs/src/{differential_forms.md => how_it_works.md} (57%) diff --git a/docs/src/differential_forms.md b/docs/src/how_it_works.md similarity index 57% rename from docs/src/differential_forms.md rename to docs/src/how_it_works.md index b8e34ab3..0a51cb9d 100644 --- a/docs/src/differential_forms.md +++ b/docs/src/how_it_works.md @@ -1,6 +1,17 @@ -# Differential Forms (Work in Progress) +# How it Works (Work in Progress) -**MeshIntegrals.jl** uses differential forms to perform numerical integration. Every supported `Meshes.Geometry` type is defined as having a parametric function that maps from a local coordinate system to every point on the geometry. For example, +This page will explain how this package works by example... + +Let $f$ be a function to be integrated throughout the volume bounded by a unit sphere. This integral is often expressed as simply +```math +\iiint f(x, y, z) ~ \text{d}V +``` + +## Parametric Functions + +**!TODO! update segment to Ball** + +Every supported `Meshes.Geometry` type is defined as having a parametric function that maps from a local coordinate system to every point on the geometry. For example, ```julia a = Meshes.Point(0, 0) b = Meshes.Point(2, 4) @@ -13,18 +24,20 @@ segment(0.5) == Point(1, 2) segment(1) == b ``` -The general solution for integrating such a geometry in 3D space can look something like the following, where $t$ is a parametric coordinate used to generate points in the domain. +... where $t$ is a parametric coordinate used to generate points in the domain. + + -**TODO: update all of the above for a 2D geometry (Sphere?) to make it more interesting/relevant. Something simple enough to follow but non-trivial.** +## Differential Forms -Using differential forms, the general solution for integrating a geometry with two parametric dimensions, $t_1$ and $t_2$, is +Using differential forms, the general solution for integrating a geometry with three parametric dimensions ($t_1$, $t_2$, and $t_3$) is ```math -\iint f(x, y, z) ~ \text{d}t_1 \wedge \text{d}t_2 +\iiint f(x, y, z) ~ \left\| \text{d}t_1 \wedge \text{d}t_2 \wedge \text{d}t_3 \right\| ``` Since `Meshes.Geometry`s parametric functions have arguments that are defined on the domain $[0,1]$, this is equivalent to ```math -\int_0^1 \int_0^1 f(x, y, z) ~ \text{d}t_1 \wedge \text{d}t_2 +\int_0^1 \int_0^1 \int_0^1 f(x, y, z) ~ \left\| \text{d}t_1 \wedge \text{d}t_2 \wedge \text{d}t_3 \right\| ``` For every point in the integration domain where the integrand function is evaluated, the differential element $\text{d}t_1 \wedge \text{d}t_2$ is calculated using [the Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of the parametric function. For a two-dimensional surface this Jacobian consists of two vectors, each pointing in the direction that the parametric function's output point will move by changing each input argument. The differential element, then, is simply the magnitude of the exterior product of these vectors. \ No newline at end of file From b00e92b390b9b9932a30f0081e946c5d305cda84 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sat, 2 Nov 2024 23:35:19 -0400 Subject: [PATCH 82/89] Add new page to directory --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index 988736cb..ad74209e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,6 +12,7 @@ makedocs( pages = [ "Home" => [ "About" => "index.md", + "How it Works" => "how_it_works.md", "Support Matrix" => "supportmatrix.md", "Example Usage" => "usage.md" ], From e4937e24521017a2509f683cdd4f16cfea98b6f4 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Sun, 3 Nov 2024 13:08:14 -0500 Subject: [PATCH 83/89] Update index page with README contents --- docs/src/index.md | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 52edc831..f947eb64 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -10,5 +10,38 @@ [![Coveralls](https://coveralls.io/repos/github/JuliaGeometry/MeshIntegrals.jl/badge.svg?branch=main)](https://coveralls.io/github/JuliaGeometry/MeshIntegrals.jl?branch=main) [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) -**MeshIntegrals.jl** is a Julia library that leverages differential forms to implement fast -and easy numerical integration of field equations over geometric domains. + +**MeshIntegrals.jl** uses differential forms to enable fast and easy numerical integration of arbitrary integrand functions over domains defined via [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) geometries. This is achieved using: +- 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...)` + +These solvers have support for integrand functions that produce scalars, vectors, and [**Unitful.jl**](https://github.com/PainterQubits/Unitful.jl) `Quantity` types. While HCubature.jl does not natively support `Quantity` type integrands, this package provides a compatibility layer to enable this feature. + +## Usage + +### Basic + +```julia +integral(f, geometry) +``` +Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry`. A default integration method will be automatically selected according to the geometry: `GaussKronrod()` for 1D, and `HAdaptiveCubature()` for all others. + +```julia +integral(f, geometry, rule) +``` +Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`. + +Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics. + +### Aliases + +```julia +lineintegral(f, geometry) +surfaceintegral(f, geometry) +volumeintegral(f, geometry) +``` +Alias functions are provided for convenience. These are simply wrappers for `integral` that also validate that the provided `geometry` has the expected number of parametric dimensions. Like with `integral`, a `rule` can also optionally be specified as a third argument. +- `lineintegral` is used for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc) +- `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc) +- `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc) From 3ae0a8b66830e92220170346bdd9659aefe846fe Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Mon, 4 Nov 2024 09:00:43 -0500 Subject: [PATCH 84/89] add EnzymeExt --- .gitignore | 4 ++++ Project.toml | 7 +++++++ ext/EnzymeExt.jl | 12 ++++++++++++ src/MeshIntegrals.jl | 4 ++-- src/differentiation.jl | 17 +++++++++++++++-- 5 files changed, 40 insertions(+), 4 deletions(-) create mode 100644 ext/EnzymeExt.jl diff --git a/.gitignore b/.gitignore index 29126e47..2335c33a 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,7 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml + +# Development related +.vscode +dev/ diff --git a/Project.toml b/Project.toml index 5a1d4711..2c8c99f5 100644 --- a/Project.toml +++ b/Project.toml @@ -12,9 +12,16 @@ Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +[weakdeps] +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + +[extensions] +EnzymeExt = "Enzyme" + [compat] CliffordNumbers = "0.1.4" CoordRefSystems = "0.12, 0.13, 0.14, 0.15" +Enzyme = "0.13.12" FastGaussQuadrature = "1" HCubature = "1.5" LinearAlgebra = "1" diff --git a/ext/EnzymeExt.jl b/ext/EnzymeExt.jl new file mode 100644 index 00000000..4af729f6 --- /dev/null +++ b/ext/EnzymeExt.jl @@ -0,0 +1,12 @@ +module EnzymeExt + +using MeshIntegrals +using Meshes: Geometry +using Enzyme + +function MeshIntegrals.jacobian( + geometry::Geometry, ts::V, ::AutoEnzyme) where {V <: Union{AbstractVector, Tuple}} + to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...)) +end + +end diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 5fbac40f..1124035b 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -1,7 +1,7 @@ module MeshIntegrals using CliffordNumbers: CliffordNumbers, VGA, ∧ using CoordRefSystems: CoordRefSystems, CRS -using Meshes: Meshes, Geometry +using Meshes: Meshes, Geometry, to import FastGaussQuadrature import HCubature @@ -10,7 +10,7 @@ import QuadGK import Unitful include("differentiation.jl") -export DifferentiationMethod, Analytical, FiniteDifference, jacobian +export DifferentiationMethod, Analytical, FiniteDifference, AutoEnzyme, jacobian include("utils.jl") diff --git a/src/differentiation.jl b/src/differentiation.jl index d0941129..d3d19ec9 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -42,8 +42,21 @@ These solutions are currently defined only for a subset of geometry types. """ struct Analytical <: DifferentiationMethod end +""" + AutoEnzyme() + +Use to specify use of the Enzyme.jl for calculating derivatives. +""" +struct AutoEnzyme <: DifferentiationMethod + function AutoEnzyme() + if :Enzyme ∉ names(Main, imported = true) + error("Please load Enzyme.jl to use AutoEnzyme().") + end + return new() + end +end + # Future Support: -# struct AutoEnzyme <: DifferentiationMethod end # struct AutoZygote <: DifferentiationMethod end ################################################################################ @@ -86,7 +99,7 @@ function jacobian( # Get the partial derivative along the n'th axis via finite difference # approximation, where ts is the current parametric position function ∂ₙr(ts, n) - # Build left/right parametric coordinates with non-allocating iterators + # Build left/right parametric coordinates with non-allocating iterators left = Iterators.map(((i, t),) -> i == n ? t - ε : t, enumerate(ts)) right = Iterators.map(((i, t),) -> i == n ? t + ε : t, enumerate(ts)) # Select orientation of finite-diff From d596ef7b51e76e93af5b023a96e7db3ccedebf1a Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Mon, 4 Nov 2024 19:02:41 -0500 Subject: [PATCH 85/89] clean up some imports and naming --- Project.toml | 2 +- ext/{EnzymeExt.jl => MeshIntegralsEnzymeExt.jl} | 8 ++++---- src/MeshIntegrals.jl | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) rename ext/{EnzymeExt.jl => MeshIntegralsEnzymeExt.jl} (60%) diff --git a/Project.toml b/Project.toml index 2c8c99f5..4e90add0 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" [extensions] -EnzymeExt = "Enzyme" +MeshIntegralsEnzymeExt = "Enzyme" [compat] CliffordNumbers = "0.1.4" diff --git a/ext/EnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl similarity index 60% rename from ext/EnzymeExt.jl rename to ext/MeshIntegralsEnzymeExt.jl index 4af729f6..bda73856 100644 --- a/ext/EnzymeExt.jl +++ b/ext/MeshIntegralsEnzymeExt.jl @@ -1,8 +1,8 @@ -module EnzymeExt +module MeshIntegralsEnzymeExt -using MeshIntegrals -using Meshes: Geometry -using Enzyme +using MeshIntegrals: MeshIntegrals, AutoEnzyme +using Meshes: Geometry, to +using Enzyme: Enzyme function MeshIntegrals.jacobian( geometry::Geometry, ts::V, ::AutoEnzyme) where {V <: Union{AbstractVector, Tuple}} diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 1124035b..66c2ec36 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -1,7 +1,7 @@ module MeshIntegrals using CliffordNumbers: CliffordNumbers, VGA, ∧ using CoordRefSystems: CoordRefSystems, CRS -using Meshes: Meshes, Geometry, to +using Meshes: Meshes, Geometry import FastGaussQuadrature import HCubature From 5a3d4ab12e64a853139c2ce537216523cc513f8d Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 11 Nov 2024 23:33:04 -0500 Subject: [PATCH 86/89] Tweak formatting and update ts type --- ext/MeshIntegralsEnzymeExt.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/ext/MeshIntegralsEnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl index bda73856..a0ad62c8 100644 --- a/ext/MeshIntegralsEnzymeExt.jl +++ b/ext/MeshIntegralsEnzymeExt.jl @@ -4,9 +4,12 @@ using MeshIntegrals: MeshIntegrals, AutoEnzyme using Meshes: Geometry, to using Enzyme: Enzyme -function MeshIntegrals.jacobian( - geometry::Geometry, ts::V, ::AutoEnzyme) where {V <: Union{AbstractVector, Tuple}} - to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...)) +function Meshes.jacobian( + geometry::Geometry, + ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}, + ::AutoEnzyme +) where {T <: AbstractFloat} + return to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...)) end end From a04ecd46696c8e726322e9995ee24b032da76ffb Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Tue, 19 Nov 2024 20:53:31 -0500 Subject: [PATCH 87/89] edit tests to use different AD backends --- ext/MeshIntegralsEnzymeExt.jl | 8 +- src/differentiation.jl | 12 +- test/Project.toml | 1 + test/combinations.jl | 462 +++++++++++++++++----------------- test/runtests.jl | 36 +++ test/utils.jl | 3 + 6 files changed, 281 insertions(+), 241 deletions(-) diff --git a/ext/MeshIntegralsEnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl index a0ad62c8..ccc6f36d 100644 --- a/ext/MeshIntegralsEnzymeExt.jl +++ b/ext/MeshIntegralsEnzymeExt.jl @@ -1,15 +1,15 @@ module MeshIntegralsEnzymeExt using MeshIntegrals: MeshIntegrals, AutoEnzyme -using Meshes: Geometry, to +using Meshes: Meshes using Enzyme: Enzyme -function Meshes.jacobian( - geometry::Geometry, +function MeshIntegrals.jacobian( + geometry::Meshes.Geometry, ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}, ::AutoEnzyme ) where {T <: AbstractFloat} - return to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...)) + return Meshes.to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...)) end end diff --git a/src/differentiation.jl b/src/differentiation.jl index d3d19ec9..53bf9071 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -48,12 +48,12 @@ struct Analytical <: DifferentiationMethod end Use to specify use of the Enzyme.jl for calculating derivatives. """ struct AutoEnzyme <: DifferentiationMethod - function AutoEnzyme() - if :Enzyme ∉ names(Main, imported = true) - error("Please load Enzyme.jl to use AutoEnzyme().") - end - return new() - end + #function AutoEnzyme() + # if :Enzyme ∉ names(Main, imported = true) + # error("Please load Enzyme.jl to use AutoEnzyme().") + # end + # return new() + #end end # Future Support: diff --git a/test/Project.toml b/test/Project.toml index 7cf2952f..9bec3654 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,7 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" diff --git a/test/combinations.jl b/test/combinations.jl index aa5e87f5..b7cdc37a 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -15,20 +15,20 @@ # Scalar integrand sol = (π - π * exp(-radius^2)) * u"m^2" - @test integral(f, ball, GaussLegendre(100)) ≈ sol - @test integral(f, ball, GaussKronrod()) ≈ sol - @test integral(f, ball, HAdaptiveCubature()) ≈ sol + @test integral_test(f, ball, GaussLegendre(100)) ≈ sol + @test integral_test(f, ball, GaussKronrod()) ≈ sol + @test integral_test(f, ball, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, ball, GaussLegendre(100)) ≈ vsol - @test integral(fv, ball, GaussKronrod()) ≈ vsol - @test integral(fv, ball, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, ball, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, ball, GaussKronrod()) ≈ vsol + @test integral_test(fv, ball, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, ball) - @test surfaceintegral(f, ball) ≈ sol - @test_throws "not supported" volumeintegral(f, ball) + @test_throws "not supported" lineintegral_test(f, ball) + @test surfaceintegral_test(f, ball) ≈ sol + @test_throws "not supported" volumeintegral_test(f, ball) end @testitem "Meshes.Ball 3D" setup=[Setup] begin @@ -49,20 +49,20 @@ end # Scalar integrand r = ustrip(u"m", radius) sol = (π^(3 / 2) * erf(r) - 2π * exp(-r^2) * r) * u"m^3" - @test integral(f, ball, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, ball, GaussKronrod())≈sol - @test integral(f, ball, HAdaptiveCubature()) ≈ sol + @test integral_test(f, ball, GaussLegendre(100)) ≈ sol + @test_throws "not supported" integral_test(f, ball, GaussKronrod())≈sol + @test integral_test(f, ball, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, ball, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, ball, GaussKronrod())≈vsol - @test integral(fv, ball, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, ball, GaussLegendre(100)) ≈ vsol + @test_throws "not supported" integral_test(fv, ball, GaussKronrod())≈vsol + @test integral_test(fv, ball, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, ball) - @test_throws "not supported" surfaceintegral(f, ball) - @test volumeintegral(f, ball) ≈ sol + @test_throws "not supported" lineintegral_test(f, ball) + @test_throws "not supported" surfaceintegral_test(f, ball) + @test volumeintegral_test(f, ball) ≈ sol end @testitem "Meshes.BezierCurve" setup=[Setup] begin @@ -78,20 +78,20 @@ end # Scalar integrand sol = 2π * u"Ω" - @test integral(f, curve, GaussLegendre(100))≈sol rtol=0.5e-2 - @test integral(f, curve, GaussKronrod())≈sol rtol=0.5e-2 - @test integral(f, curve, HAdaptiveCubature())≈sol rtol=0.5e-2 + @test integral_test(f, curve, GaussLegendre(100))≈sol rtol=0.5e-2 + @test integral_test(f, curve, GaussKronrod())≈sol rtol=0.5e-2 + @test integral_test(f, curve, HAdaptiveCubature())≈sol rtol=0.5e-2 # Vector integrand vsol = fill(sol, 3) - @test integral(fv, curve, GaussLegendre(100))≈vsol rtol=0.5e-2 - @test integral(fv, curve, GaussKronrod())≈vsol rtol=0.5e-2 - @test integral(fv, curve, HAdaptiveCubature())≈vsol rtol=0.5e-2 + @test integral_test(fv, curve, GaussLegendre(100))≈vsol rtol=0.5e-2 + @test integral_test(fv, curve, GaussKronrod())≈vsol rtol=0.5e-2 + @test integral_test(fv, curve, HAdaptiveCubature())≈vsol rtol=0.5e-2 # Integral aliases - @test lineintegral(f, curve)≈sol rtol=0.5e-2 - @test_throws "not supported" surfaceintegral(f, curve) - @test_throws "not supported" volumeintegral(f, curve) + @test lineintegral_test(f, curve)≈sol rtol=0.5e-2 + @test_throws "not supported" surfaceintegral_test(f, curve) + @test_throws "not supported" volumeintegral_test(f, curve) # Check Bezier-specific jacobian bounds @test_throws DomainError jacobian(curve, [1.1]) @@ -109,20 +109,20 @@ end # Scalar integrand sol = π * a^2 / 4 * u"Ω" - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(f, box, GaussKronrod()) ≈ sol - @test integral(f, box, HAdaptiveCubature()) ≈ sol + @test integral_test(f, box, GaussLegendre(100))≈sol rtol=1e-6 + @test integral_test(f, box, GaussKronrod()) ≈ sol + @test integral_test(f, box, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, box, GaussKronrod()) ≈ vsol - @test integral(fv, box, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 + @test integral_test(fv, box, GaussKronrod()) ≈ vsol + @test integral_test(fv, box, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, box) ≈ sol - @test_throws "not supported" surfaceintegral(f, box) - @test_throws "not supported" volumeintegral(f, box) + @test lineintegral_test(f, box) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, box) + @test_throws "not supported" volumeintegral_test(f, box) end @testitem "Meshes.Box 2D" setup=[Setup] begin @@ -137,20 +137,20 @@ end # Scalar integrand sol = 2a * (π * a^2 / 4) * u"Ω" - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(f, box, GaussKronrod()) ≈ sol - @test integral(f, box, HAdaptiveCubature()) ≈ sol + @test integral_test(f, box, GaussLegendre(100))≈sol rtol=1e-6 + @test integral_test(f, box, GaussKronrod()) ≈ sol + @test integral_test(f, box, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, box, GaussKronrod()) ≈ vsol - @test integral(fv, box, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 + @test integral_test(fv, box, GaussKronrod()) ≈ vsol + @test integral_test(fv, box, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test surfaceintegral(f, box) ≈ sol - @test_throws "not supported" volumeintegral(f, box) + @test_throws "not supported" lineintegral_test(f, box) + @test surfaceintegral_test(f, box) ≈ sol + @test_throws "not supported" volumeintegral_test(f, box) # Test jacobian with wrong number of parametric coordinates @test_throws ArgumentError jacobian(box, zeros(3)) @@ -168,20 +168,20 @@ end # Scalar integrand sol = 3a^2 * (π * a^2 / 4) * u"Ω" - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test_throws "not supported" integral(f, box, GaussKronrod()) - @test integral(f, box, HAdaptiveCubature()) ≈ sol + @test integral_test(f, box, GaussLegendre(100))≈sol rtol=1e-6 + @test_throws "not supported" integral_test(f, box, GaussKronrod()) + @test integral_test(f, box, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test_throws "not supported" integral(fv, box, GaussKronrod()) - @test integral(fv, box, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 + @test_throws "not supported" integral_test(fv, box, GaussKronrod()) + @test integral_test(fv, box, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test_throws "not supported" surfaceintegral(f, box) - @test volumeintegral(f, box) ≈ sol + @test_throws "not supported" lineintegral_test(f, box) + @test_throws "not supported" surfaceintegral_test(f, box) + @test volumeintegral_test(f, box) ≈ sol end @testitem "Meshes.Box 4D" tags=[:extended] setup=[Setup] begin @@ -197,20 +197,20 @@ end # Scalar integrand sol = 4a^3 * (π * a^2 / 4) * u"Ω" - @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6 - @test_throws "not supported" integral(f, box, GaussKronrod()) - @test integral(f, box, HAdaptiveCubature(rtol = 1e-6))≈sol rtol=1e-6 + @test integral_test(f, box, GaussLegendre(100))≈sol rtol=1e-6 + @test_throws "not supported" integral_test(f, box, GaussKronrod()) + @test integral_test(f, box, HAdaptiveCubature(rtol = 1e-6))≈sol rtol=1e-6 # Vector integrand vsol = fill(sol, 3) - @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 - @test_throws "not supported" integral(fv, box, GaussKronrod()) - @test integral(fv, box, HAdaptiveCubature(rtol = 1e-6))≈vsol rtol=1e-6 + @test integral_test(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 + @test_throws "not supported" integral_test(fv, box, GaussKronrod()) + @test integral_test(fv, box, HAdaptiveCubature(rtol = 1e-6))≈vsol rtol=1e-6 # Integral aliases - @test_throws "not supported" lineintegral(f, box) - @test_throws "not supported" surfaceintegral(f, box) - @test_throws "not supported" volumeintegral(f, box) + @test_throws "not supported" lineintegral_test(f, box) + @test_throws "not supported" surfaceintegral_test(f, box) + @test_throws "not supported" volumeintegral_test(f, box) end @testitem "Meshes.Circle" setup=[Setup] begin @@ -229,20 +229,20 @@ end # Scalar integrand sol = 2π * radius * exp(-radius^2) * u"m" - @test integral(f, circle, GaussLegendre(100)) ≈ sol - @test integral(f, circle, GaussKronrod()) ≈ sol - @test integral(f, circle, HAdaptiveCubature()) ≈ sol + @test integral_test(f, circle, GaussLegendre(100)) ≈ sol + @test integral_test(f, circle, GaussKronrod()) ≈ sol + @test integral_test(f, circle, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, circle, GaussLegendre(100)) ≈ vsol - @test integral(fv, circle, GaussKronrod()) ≈ vsol - @test integral(fv, circle, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, circle, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, circle, GaussKronrod()) ≈ vsol + @test integral_test(fv, circle, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, circle) ≈ sol - @test_throws "not supported" surfaceintegral(f, circle) - @test_throws "not supported" volumeintegral(f, circle) + @test lineintegral_test(f, circle) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, circle) + @test_throws "not supported" volumeintegral_test(f, circle) end @testitem "Meshes.Cone" setup=[Setup] begin @@ -261,20 +261,20 @@ end # Scalar integrand sol = _volume_cone_rightcircular(r, h) - @test integral(f, cone, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, cone, GaussKronrod()) - @test integral(f, cone, HAdaptiveCubature()) ≈ sol + @test integral_test(f, cone, GaussLegendre(100)) ≈ sol + @test_throws "not supported" integral_test(f, cone, GaussKronrod()) + @test integral_test(f, cone, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, cone, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, cone, GaussKronrod()) - @test integral(fv, cone, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, cone, GaussLegendre(100)) ≈ vsol + @test_throws "not supported" integral_test(fv, cone, GaussKronrod()) + @test integral_test(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 + @test_throws "not supported" lineintegral_test(f, cone) + @test_throws "not supported" surfaceintegral_test(f, cone) + @test volumeintegral_test(f, cone) ≈ sol end @testitem "Meshes.ConeSurface" setup=[Setup] begin @@ -293,20 +293,20 @@ end # Scalar integrand sol = _area_cone_rightcircular(h, r) - @test integral(f, cone, GaussLegendre(100))≈sol rtol=1e-6 - @test integral(f, cone, GaussKronrod())≈sol rtol=1e-6 - @test integral(f, cone, HAdaptiveCubature()) ≈ sol + @test integral_test(f, cone, GaussLegendre(100))≈sol rtol=1e-6 + @test integral_test(f, cone, GaussKronrod())≈sol rtol=1e-6 + @test integral_test(f, cone, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, cone, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, cone, GaussKronrod())≈vsol rtol=1e-6 - @test integral(fv, cone, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, cone, GaussLegendre(100))≈vsol rtol=1e-6 + @test integral_test(fv, cone, GaussKronrod())≈vsol rtol=1e-6 + @test integral_test(fv, cone, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, cone) - @test surfaceintegral(f, cone) ≈ sol - @test_throws "not supported" volumeintegral(f, cone) + @test_throws "not supported" lineintegral_test(f, cone) + @test surfaceintegral_test(f, cone) ≈ sol + @test_throws "not supported" volumeintegral_test(f, cone) end @testitem "Meshes.Cylinder" setup=[Setup] begin @@ -319,20 +319,20 @@ end # Scalar integrand sol = Meshes.measure(cyl) - @test integral(f, cyl, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, cyl, GaussKronrod()) - @test integral(f, cyl, HAdaptiveCubature()) ≈ sol + @test integral_test(f, cyl, GaussLegendre(100)) ≈ sol + @test_throws "not supported" integral_test(f, cyl, GaussKronrod()) + @test integral_test(f, cyl, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, cyl, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, cyl, GaussKronrod()) - @test integral(fv, cyl, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, cyl, GaussLegendre(100)) ≈ vsol + @test_throws "not supported" integral_test(fv, cyl, GaussKronrod()) + @test integral_test(fv, cyl, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, cyl) - @test_throws "not supported" surfaceintegral(f, cyl) - @test volumeintegral(f, cyl) ≈ sol + @test_throws "not supported" lineintegral_test(f, cyl) + @test_throws "not supported" surfaceintegral_test(f, cyl) + @test volumeintegral_test(f, cyl) ≈ sol end @testitem "Meshes.CylinderSurface" setup=[Setup] begin @@ -345,20 +345,20 @@ end # Scalar integrand sol = Meshes.measure(cyl) - @test integral(f, cyl, GaussLegendre(100)) ≈ sol - @test integral(f, cyl, GaussKronrod()) ≈ sol - @test integral(f, cyl, HAdaptiveCubature()) ≈ sol + @test integral_test(f, cyl, GaussLegendre(100)) ≈ sol + @test integral_test(f, cyl, GaussKronrod()) ≈ sol + @test integral_test(f, cyl, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, cyl, GaussLegendre(100)) ≈ vsol - @test integral(fv, cyl, GaussKronrod()) ≈ vsol - @test integral(fv, cyl, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, cyl, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, cyl, GaussKronrod()) ≈ vsol + @test integral_test(fv, cyl, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, cyl) - @test surfaceintegral(f, cyl) ≈ sol - @test_throws "not supported" volumeintegral(f, cyl) + @test_throws "not supported" lineintegral_test(f, cyl) + @test surfaceintegral_test(f, cyl) ≈ sol + @test_throws "not supported" volumeintegral_test(f, cyl) end @testitem "Meshes.Disk" setup=[Setup] begin @@ -377,20 +377,20 @@ end # Scalar integrand sol = (π - π * exp(-radius^2)) * u"m^2" - @test integral(f, disk, GaussLegendre(100)) ≈ sol - @test integral(f, disk, GaussKronrod()) ≈ sol - @test integral(f, disk, HAdaptiveCubature()) ≈ sol + @test integral_test(f, disk, GaussLegendre(100)) ≈ sol + @test integral_test(f, disk, GaussKronrod()) ≈ sol + @test integral_test(f, disk, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, disk, GaussLegendre(100)) ≈ vsol - @test integral(fv, disk, GaussKronrod()) ≈ vsol - @test integral(fv, disk, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, disk, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, disk, GaussKronrod()) ≈ vsol + @test integral_test(fv, disk, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, disk) - @test surfaceintegral(f, disk) ≈ sol - @test_throws "not supported" volumeintegral(f, disk) + @test_throws "not supported" lineintegral_test(f, disk) + @test surfaceintegral_test(f, disk) ≈ sol + @test_throws "not supported" volumeintegral_test(f, disk) end @testitem "Meshes.Ellipsoid" setup=[Setup] begin @@ -404,20 +404,20 @@ end # Tolerances are higher due to `measure` being only an approximation # Scalar integrand sol = Meshes.measure(ellipsoid) - @test integral(f, ellipsoid, GaussLegendre(100))≈sol rtol=1e-2 - @test integral(f, ellipsoid, GaussKronrod())≈sol rtol=1e-2 - @test integral(f, ellipsoid, HAdaptiveCubature())≈sol rtol=1e-2 + @test integral_test(f, ellipsoid, GaussLegendre(100))≈sol rtol=1e-2 + @test integral_test(f, ellipsoid, GaussKronrod())≈sol rtol=1e-2 + @test integral_test(f, ellipsoid, HAdaptiveCubature())≈sol rtol=1e-2 # Vector integrand vsol = fill(sol, 3) - @test integral(fv, ellipsoid, GaussLegendre(100))≈vsol rtol=1e-2 - @test integral(fv, ellipsoid, GaussKronrod())≈vsol rtol=1e-2 - @test integral(fv, ellipsoid, HAdaptiveCubature())≈vsol rtol=1e-2 + @test integral_test(fv, ellipsoid, GaussLegendre(100))≈vsol rtol=1e-2 + @test integral_test(fv, ellipsoid, GaussKronrod())≈vsol rtol=1e-2 + @test integral_test(fv, ellipsoid, HAdaptiveCubature())≈vsol rtol=1e-2 # Integral aliases - @test_throws "not supported" lineintegral(f, ellipsoid) - @test surfaceintegral(f, ellipsoid)≈sol rtol=1e-2 - @test_throws "not supported" volumeintegral(f, ellipsoid) + @test_throws "not supported" lineintegral_test(f, ellipsoid) + @test surfaceintegral_test(f, ellipsoid)≈sol rtol=1e-2 + @test_throws "not supported" volumeintegral_test(f, ellipsoid) end @testitem "Meshes.FrustumSurface" setup=[Setup] begin @@ -448,15 +448,15 @@ 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 rtol=1e-6 - @test integral(f, frustum, GaussKronrod())≈sol rtol=1e-6 - @test integral(f, frustum, HAdaptiveCubature()) ≈ sol + @test integral_test(f, frustum, GaussLegendre(100))≈sol rtol=1e-6 + @test integral_test(f, frustum, GaussKronrod())≈sol rtol=1e-6 + @test integral_test(f, frustum, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, frustum, GaussLegendre(100))≈vsol rtol=1e-6 - @test integral(fv, frustum, GaussKronrod())≈vsol rtol=1e-6 - @test integral(fv, frustum, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, frustum, GaussLegendre(100))≈vsol rtol=1e-6 + @test integral_test(fv, frustum, GaussKronrod())≈vsol rtol=1e-6 + @test integral_test(fv, frustum, HAdaptiveCubature()) ≈ vsol end @testitem "Meshes.Hexahedron" setup=[Setup] begin @@ -468,20 +468,20 @@ end # Scalar integrand sol = Meshes.measure(hexahedron) - @test integral(f, hexahedron, GaussLegendre(100)) ≈ sol - @test_throws "not supported" integral(f, hexahedron, GaussKronrod())≈sol - @test integral(f, hexahedron, HAdaptiveCubature()) ≈ sol + @test integral_test(f, hexahedron, GaussLegendre(100)) ≈ sol + @test_throws "not supported" integral_test(f, hexahedron, GaussKronrod())≈sol + @test integral_test(f, hexahedron, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, hexahedron, GaussLegendre(100)) ≈ vsol - @test_throws "not supported" integral(fv, hexahedron, GaussKronrod())≈vsol - @test integral(fv, hexahedron, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, hexahedron, GaussLegendre(100)) ≈ vsol + @test_throws "not supported" integral_test(fv, hexahedron, GaussKronrod())≈vsol + @test integral_test(fv, hexahedron, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, hexahedron) - @test_throws "not supported" surfaceintegral(f, hexahedron) - @test volumeintegral(f, hexahedron) ≈ sol + @test_throws "not supported" lineintegral_test(f, hexahedron) + @test_throws "not supported" surfaceintegral_test(f, hexahedron) + @test volumeintegral_test(f, hexahedron) ≈ sol end @testitem "Meshes.Line" setup=[Setup] begin @@ -522,20 +522,20 @@ end # Scalar integrand sol = Meshes.measure(parab) - @test integral(f, parab, GaussLegendre(100)) ≈ sol - @test integral(f, parab, GaussKronrod()) ≈ sol - @test integral(f, parab, HAdaptiveCubature()) ≈ sol + @test integral_test(f, parab, GaussLegendre(100)) ≈ sol + @test integral_test(f, parab, GaussKronrod()) ≈ sol + @test integral_test(f, parab, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, parab, GaussLegendre(100)) ≈ vsol - @test integral(fv, parab, GaussKronrod()) ≈ vsol - @test integral(fv, parab, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, parab, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, parab, GaussKronrod()) ≈ vsol + @test integral_test(fv, parab, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, parab) - @test surfaceintegral(f, parab) ≈ sol - @test_throws "not supported" volumeintegral(f, parab) + @test_throws "not supported" lineintegral_test(f, parab) + @test surfaceintegral_test(f, parab) ≈ sol + @test_throws "not supported" volumeintegral_test(f, parab) end @testitem "ParametrizedCurve" setup=[Setup] begin @@ -559,29 +559,29 @@ end # Scalar integrand sol = 2π * radius * exp(-radius^2) * u"m" - @test integral(f, curve_cart, GaussLegendre(100)) ≈ sol - @test integral(f, curve_cart, GaussKronrod()) ≈ sol - @test integral(f, curve_cart, HAdaptiveCubature()) ≈ sol - @test integral(f, curve_polar, GaussLegendre(100)) ≈ sol - @test integral(f, curve_polar, GaussKronrod()) ≈ sol - @test integral(f, curve_polar, HAdaptiveCubature()) ≈ sol + @test integral_test(f, curve_cart, GaussLegendre(100)) ≈ sol + @test integral_test(f, curve_cart, GaussKronrod()) ≈ sol + @test integral_test(f, curve_cart, HAdaptiveCubature()) ≈ sol + @test integral_test(f, curve_polar, GaussLegendre(100)) ≈ sol + @test integral_test(f, curve_polar, GaussKronrod()) ≈ sol + @test integral_test(f, curve_polar, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, curve_cart, GaussLegendre(100)) ≈ vsol - @test integral(fv, curve_cart, GaussKronrod()) ≈ vsol - @test integral(fv, curve_cart, HAdaptiveCubature()) ≈ vsol - @test integral(fv, curve_polar, GaussLegendre(100)) ≈ vsol - @test integral(fv, curve_polar, GaussKronrod()) ≈ vsol - @test integral(fv, curve_polar, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, curve_cart, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, curve_cart, GaussKronrod()) ≈ vsol + @test integral_test(fv, curve_cart, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, curve_polar, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, curve_polar, GaussKronrod()) ≈ vsol + @test integral_test(fv, curve_polar, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, curve_cart) ≈ sol - @test_throws "not supported" surfaceintegral(f, curve_cart) - @test_throws "not supported" volumeintegral(f, curve_cart) - @test lineintegral(f, curve_polar) ≈ sol - @test_throws "not supported" surfaceintegral(f, curve_polar) - @test_throws "not supported" volumeintegral(f, curve_polar) + @test lineintegral_test(f, curve_cart) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, curve_cart) + @test_throws "not supported" volumeintegral_test(f, curve_cart) + @test lineintegral_test(f, curve_polar) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, curve_polar) + @test_throws "not supported" volumeintegral_test(f, curve_polar) end end @@ -626,20 +626,20 @@ end # Scalar integrand sol = 0.5 * π * erf(1)^2 * u"m^2" - @test integral(f, quadrangle, GaussLegendre(100)) ≈ sol - @test integral(f, quadrangle, GaussKronrod()) ≈ sol - @test integral(f, quadrangle, HAdaptiveCubature()) ≈ sol + @test integral_test(f, quadrangle, GaussLegendre(100)) ≈ sol + @test integral_test(f, quadrangle, GaussKronrod()) ≈ sol + @test integral_test(f, quadrangle, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, quadrangle, GaussLegendre(100)) ≈ vsol - @test integral(fv, quadrangle, GaussKronrod()) ≈ vsol - @test integral(fv, quadrangle, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, quadrangle, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, quadrangle, GaussKronrod()) ≈ vsol + @test integral_test(fv, quadrangle, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, quadrangle) - @test surfaceintegral(f, quadrangle) ≈ sol - @test_throws "not supported" volumeintegral(f, quadrangle) + @test_throws "not supported" lineintegral_test(f, quadrangle) + @test surfaceintegral_test(f, quadrangle) ≈ sol + @test_throws "not supported" volumeintegral_test(f, quadrangle) end @testitem "Meshes.Ray" setup=[Setup] begin @@ -686,20 +686,20 @@ end # Scalar integrand sol = 14.0u"A" - @test integral(f, rope, GaussLegendre(100)) ≈ sol - @test integral(f, rope, GaussKronrod()) ≈ sol - @test integral(f, rope, HAdaptiveCubature()) ≈ sol + @test integral_test(f, rope, GaussLegendre(100)) ≈ sol + @test integral_test(f, rope, GaussKronrod()) ≈ sol + @test integral_test(f, rope, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, rope, GaussLegendre(100)) ≈ vsol - @test integral(fv, rope, GaussKronrod()) ≈ vsol - @test integral(fv, rope, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, rope, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, rope, GaussKronrod()) ≈ vsol + @test integral_test(fv, rope, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, rope) ≈ sol - @test_throws "not supported" surfaceintegral(f, rope) - @test_throws "not supported" volumeintegral(f, rope) + @test lineintegral_test(f, rope) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, rope) + @test_throws "not supported" volumeintegral_test(f, rope) end @testitem "Meshes.Rope" setup=[Setup] begin @@ -717,20 +717,20 @@ end # Scalar integrand sol = 7.0u"A" - @test integral(f, rope, GaussLegendre(100)) ≈ sol - @test integral(f, rope, GaussKronrod()) ≈ sol - @test integral(f, rope, HAdaptiveCubature()) ≈ sol + @test integral_test(f, rope, GaussLegendre(100)) ≈ sol + @test integral_test(f, rope, GaussKronrod()) ≈ sol + @test integral_test(f, rope, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, rope, GaussLegendre(100)) ≈ vsol - @test integral(fv, rope, GaussKronrod()) ≈ vsol - @test integral(fv, rope, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, rope, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, rope, GaussKronrod()) ≈ vsol + @test integral_test(fv, rope, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, rope) ≈ sol - @test_throws "not supported" surfaceintegral(f, rope) - @test_throws "not supported" volumeintegral(f, rope) + @test lineintegral_test(f, rope) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, rope) + @test_throws "not supported" volumeintegral_test(f, rope) end @testitem "Meshes.Segment" setup=[Setup] begin @@ -750,20 +750,20 @@ end # Scalar integrand sol = (a - b) / (log(a) - log(b)) * u"m" - @test integral(f, segment, GaussLegendre(100)) ≈ sol - @test integral(f, segment, GaussKronrod()) ≈ sol - @test integral(f, segment, HAdaptiveCubature()) ≈ sol + @test integral_test(f, segment, GaussLegendre(100)) ≈ sol + @test integral_test(f, segment, GaussKronrod()) ≈ sol + @test integral_test(f, segment, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, segment, GaussLegendre(100)) ≈ vsol - @test integral(fv, segment, GaussKronrod()) ≈ vsol - @test integral(fv, segment, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, segment, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, segment, GaussKronrod()) ≈ vsol + @test integral_test(fv, segment, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, segment) ≈ sol - @test_throws "not supported" surfaceintegral(f, segment) - @test_throws "not supported" volumeintegral(f, segment) + @test lineintegral_test(f, segment) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, segment) + @test_throws "not supported" volumeintegral_test(f, segment) end @testitem "Meshes.Sphere 2D" setup=[Setup] begin @@ -779,20 +779,20 @@ end # Scalar integrand sol = 2π * radius * exp(-radius^2) * u"m" - @test integral(f, sphere, GaussLegendre(100)) ≈ sol - @test integral(f, sphere, GaussKronrod()) ≈ sol - @test integral(f, sphere, HAdaptiveCubature()) ≈ sol + @test integral_test(f, sphere, GaussLegendre(100)) ≈ sol + @test integral_test(f, sphere, GaussKronrod()) ≈ sol + @test integral_test(f, sphere, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, sphere, GaussLegendre(100)) ≈ vsol - @test integral(fv, sphere, GaussKronrod()) ≈ vsol - @test integral(fv, sphere, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, sphere, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, sphere, GaussKronrod()) ≈ vsol + @test integral_test(fv, sphere, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test lineintegral(f, sphere) ≈ sol - @test_throws "not supported" surfaceintegral(f, sphere) - @test_throws "not supported" volumeintegral(f, sphere) + @test lineintegral_test(f, sphere) ≈ sol + @test_throws "not supported" surfaceintegral_test(f, sphere) + @test_throws "not supported" volumeintegral_test(f, sphere) end @testitem "Meshes.Sphere 3D" setup=[Setup] begin @@ -813,20 +813,20 @@ end # Scalar integrand sol = (2π * radius^2) + (4π / 3 * radius^2) - @test integral(f, sphere, GaussLegendre(100)) ≈ sol - @test integral(f, sphere, GaussKronrod()) ≈ sol - @test integral(f, sphere, HAdaptiveCubature()) ≈ sol + @test integral_test(f, sphere, GaussLegendre(100)) ≈ sol + @test integral_test(f, sphere, GaussKronrod()) ≈ sol + @test integral_test(f, sphere, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, sphere, GaussLegendre(100)) ≈ vsol - @test integral(fv, sphere, GaussKronrod()) ≈ vsol - @test integral(fv, sphere, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, sphere, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, sphere, GaussKronrod()) ≈ vsol + @test integral_test(fv, sphere, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, sphere) - @test surfaceintegral(f, sphere) ≈ sol - @test_throws "not supported" volumeintegral(f, sphere) + @test_throws "not supported" lineintegral_test(f, sphere) + @test surfaceintegral_test(f, sphere) ≈ sol + @test_throws "not supported" volumeintegral_test(f, sphere) end @testitem "Meshes.Tetrahedron" tags=[:extended] setup=[Setup] begin @@ -867,20 +867,20 @@ end # Scalar integrand sol = Meshes.measure(torus) - @test integral(f, torus, GaussLegendre(100)) ≈ sol - @test integral(f, torus, GaussKronrod()) ≈ sol - @test integral(f, torus, HAdaptiveCubature()) ≈ sol + @test integral_test(f, torus, GaussLegendre(100)) ≈ sol + @test integral_test(f, torus, GaussKronrod()) ≈ sol + @test integral_test(f, torus, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test integral(fv, torus, GaussLegendre(100)) ≈ vsol - @test integral(fv, torus, GaussKronrod()) ≈ vsol - @test integral(fv, torus, HAdaptiveCubature()) ≈ vsol + @test integral_test(fv, torus, GaussLegendre(100)) ≈ vsol + @test integral_test(fv, torus, GaussKronrod()) ≈ vsol + @test integral_test(fv, torus, HAdaptiveCubature()) ≈ vsol # Integral aliases - @test_throws "not supported" lineintegral(f, torus) - @test surfaceintegral(f, torus) ≈ sol - @test_throws "not supported" volumeintegral(f, torus) + @test_throws "not supported" lineintegral_test(f, torus) + @test surfaceintegral_test(f, torus) ≈ sol + @test_throws "not supported" volumeintegral_test(f, torus) end @testitem "Meshes.Triangle" setup=[Setup] begin diff --git a/test/runtests.jl b/test/runtests.jl index 380bff2d..22bcd45b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,4 +8,40 @@ using TestItems using LinearAlgebra: norm using Meshes using Unitful + + # environment settings + isCI = "CI" ∈ keys(ENV) + + # float settings + AD = if isCI + if ENV["AD"] == "FiniteDifference" + FiniteDifference() + elseif ENV["AD"] == "Enzyme" + using Enzyme: Enzyme + AutoEnzyme() + end + else + using Enzyme: Enzyme + AutoEnzyme() + end + + @info "Running tests with autodiff backend $AD" + + for g in (:integral, :lineintegral, :surfaceintegral, :volumeintegral) + newname = Symbol(g, "_test") + @eval begin + function $newname(f, geometry, rule; kwargs...) + MeshIntegrals.$g(f, geometry, rule; diff_method = AD, kwargs...) + end + end + end + + for g in (:lineintegral, :surfaceintegral, :volumeintegral) + newname = Symbol(g, "_test") + @eval begin + function $newname(f, geometry; kwargs...) + MeshIntegrals.$g(f, geometry; diff_method = AD, kwargs...) + end + end + end end diff --git a/test/utils.jl b/test/utils.jl index 2a2d8a68..dd1713e8 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -40,4 +40,7 @@ end # FiniteDifference @test FiniteDifference().ε ≈ 1e-6 + + # AutoEnzyme + @test AutoEnzyme() isa AutoEnzyme end From 87028a6b665093d0fd8825f2a0e7e4c3e9ecd18f Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Tue, 19 Nov 2024 21:16:07 -0500 Subject: [PATCH 88/89] remove nice user error message if Enzyme isn't loaded. There is some bug. --- src/differentiation.jl | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 53bf9071..3d15787b 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -47,14 +47,7 @@ struct Analytical <: DifferentiationMethod end Use to specify use of the Enzyme.jl for calculating derivatives. """ -struct AutoEnzyme <: DifferentiationMethod - #function AutoEnzyme() - # if :Enzyme ∉ names(Main, imported = true) - # error("Please load Enzyme.jl to use AutoEnzyme().") - # end - # return new() - #end -end +struct AutoEnzyme <: DifferentiationMethod end # Future Support: # struct AutoZygote <: DifferentiationMethod end From 8c2cd5e89084c0b79902c6e6b3c12a5f9a931ebc Mon Sep 17 00:00:00 2001 From: Kyle Beggs Date: Tue, 19 Nov 2024 21:17:45 -0500 Subject: [PATCH 89/89] combine test integral function wrappers in Setup --- test/runtests.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 22bcd45b..921c8d13 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,12 +33,6 @@ using TestItems function $newname(f, geometry, rule; kwargs...) MeshIntegrals.$g(f, geometry, rule; diff_method = AD, kwargs...) end - end - end - - for g in (:lineintegral, :surfaceintegral, :volumeintegral) - newname = Symbol(g, "_test") - @eval begin function $newname(f, geometry; kwargs...) MeshIntegrals.$g(f, geometry; diff_method = AD, kwargs...) end