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 7e9cfe3b..4e90add0 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" @@ -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] +MeshIntegralsEnzymeExt = "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/README.md b/README.md index 68cc0ce0..1e51d48c 100644 --- a/README.md +++ b/README.md @@ -20,27 +20,31 @@ These solvers have support for integrand functions that produce scalars, vectors ## 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, 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. +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 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) +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 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 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" ], diff --git a/docs/src/how_it_works.md b/docs/src/how_it_works.md new file mode 100644 index 00000000..0a51cb9d --- /dev/null +++ b/docs/src/how_it_works.md @@ -0,0 +1,43 @@ +# How it Works (Work in Progress) + +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) +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 +``` + +... where $t$ is a parametric coordinate used to generate points in the domain. + + + +## Differential Forms + +Using differential forms, the general solution for integrating a geometry with three parametric dimensions ($t_1$, $t_2$, and $t_3$) is +```math +\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 \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 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) diff --git a/ext/MeshIntegralsEnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl new file mode 100644 index 00000000..ccc6f36d --- /dev/null +++ b/ext/MeshIntegralsEnzymeExt.jl @@ -0,0 +1,15 @@ +module MeshIntegralsEnzymeExt + +using MeshIntegrals: MeshIntegrals, AutoEnzyme +using Meshes: Meshes +using Enzyme: Enzyme + +function MeshIntegrals.jacobian( + geometry::Meshes.Geometry, + ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}, + ::AutoEnzyme +) where {T <: AbstractFloat} + return Meshes.to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...)) +end + +end diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 19228769..66c2ec36 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -9,10 +9,10 @@ import LinearAlgebra import QuadGK import Unitful -include("utils.jl") - include("differentiation.jl") -export jacobian +export DifferentiationMethod, Analytical, FiniteDifference, AutoEnzyme, jacobian + +include("utils.jl") include("integration_rules.jl") export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature diff --git a/src/differentiation.jl b/src/differentiation.jl index acb580f1..3d15787b 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -1,22 +1,85 @@ ################################################################################ -# Jacobian and Differential Elements +# DifferentiationMethods ################################################################################ """ - jacobian(geometry, ts; ε=1e-6) + DifferentiationMethod -Calculate the Jacobian of a geometry at some parametric point `ts` using a simple -finite-difference approximation with step size `ε`. +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 + +""" + FiniteDifference(ε=1e-6) + +Use to specify use of a finite-difference approximation method with a step size +of `ε` for calculating derivatives. +""" +struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod + ε::T +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 + +""" + AutoEnzyme() + +Use to specify use of the Enzyme.jl for calculating derivatives. +""" +struct AutoEnzyme <: DifferentiationMethod end + +# Future Support: +# struct AutoZygote <: DifferentiationMethod end + +################################################################################ +# Jacobian +################################################################################ + +""" + 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 `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 -- `ε`: step size to use for the finite-difference approximation +- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use """ +function jacobian( + geometry::G, + ts::V +) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} + return jacobian(geometry, ts, _default_method(G)) +end + function jacobian( geometry::Geometry, - ts::V; - ε = 1e-6 + ts::V, + diff_method::FiniteDifference ) where {V <: Union{AbstractVector, Tuple}} Dim = Meshes.paramdim(geometry) if Dim != length(ts) @@ -24,12 +87,12 @@ function jacobian( end T = eltype(ts) - ε = T(ε) + ε = T(diff_method.ε) # 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 @@ -48,52 +111,30 @@ 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 ################################################################################ """ - differential(geometry, ts) + 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 `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`](@ref) to use """ function differential( - geometry::Geometry, - ts::V -) where {V <: Union{AbstractVector, Tuple}} + geometry::G, + ts::V, + diff_method::DifferentiationMethod = _default_method(G) +) where {G <: Geometry, V <: Union{AbstractVector, Tuple}} # Calculate the Jacobian, convert Vec -> KVector - J = jacobian(geometry, ts) + 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 d4925db7..f0f9389a 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - integral(f, geometry[, rule]; 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 @@ -12,12 +12,13 @@ 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) -- `FP`: optionally, the floating point precision desired (`Float64` by default) +- `rule`: optionally, the `IntegrationRule` used for integration (by default +`GaussKronrod()` in 1D and `HAdaptiveCubature()` else) -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 @@ -59,8 +60,9 @@ function _integral( f, geometry, rule::GaussLegendre; + diff_method::DM = _default_method(geometry), 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 +75,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, diff_method) end return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes)) @@ -84,11 +86,12 @@ function _integral( f, geometry, rule::HAdaptiveCubature; + diff_method::DM = _default_method(geometry), FP::Type{T} = Float64 -) where {T <: AbstractFloat} +) where {DM <: DifferentiationMethod, T <: AbstractFloat} N = Meshes.paramdim(geometry) - integrand(ts) = f(geometry(ts...)) * differential(geometry, ts) + 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 @@ -111,9 +114,10 @@ function _integral_gk_1d( f, geometry, rule::GaussKronrod; + diff_method::DM = _default_method(geometry), FP::Type{T} = Float64 -) where {T <: AbstractFloat} - integrand(t) = f(geometry(t)) * differential(geometry, (t,)) +) where {DM <: DifferentiationMethod, T <: AbstractFloat} + integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end @@ -121,9 +125,10 @@ function _integral_gk_2d( f, geometry2d, rule::GaussKronrod; + diff_method::DM = _default_method(geometry2d), FP::Type{T} = Float64 -) where {T <: AbstractFloat} - integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v)) +) where {DM <: DifferentiationMethod, T <: AbstractFloat} + 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 c40b3366..a6b06828 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -9,31 +9,43 @@ # keyword argument and pass through the specified algorithm choice. ################################################################################ +################################################################################ +# integral +################################################################################ """ integral(f, curve::BezierCurve, rule = GaussKronrod(); - 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 by `curve`. + +# Arguments +- `f`: an integrand function with a method `f(::Meshes.Point)` +- `curve`: a `BezierCurve` that defines the integration domain +- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration + +# 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 = _default_method(curve), 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),), 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)) @@ -43,11 +55,12 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::GaussKronrod; + diff_method::DM = _default_method(curve), 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,), diff_method) return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] end @@ -55,11 +68,12 @@ function integral( f::F, curve::Meshes.BezierCurve, rule::HAdaptiveCubature; + diff_method::DM = _default_method(curve), 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, diff_method) # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f @@ -73,3 +87,40 @@ function integral( # Reapply units return value .* integrandunits end + +################################################################################ +# jacobian +################################################################################ + +function jacobian( + 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 + 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 + +_has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index fbb398c9..655278e1 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -11,8 +11,11 @@ function integral( f::F, line::Meshes.Line, rule::GaussLegendre; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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) @@ -24,8 +27,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 + 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 @@ -33,14 +36,17 @@ function integral( f::F, line::Meshes.Line, rule::GaussKronrod; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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)) # 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 +54,11 @@ function integral( f::F, line::Meshes.Line, rule::HAdaptiveCubature; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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)) @@ -58,8 +67,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 + 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 @@ -68,8 +77,14 @@ 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 end + +################################################################################ +# jacobian +################################################################################ + +_has_analytical(::Type{T}) where {T <: Meshes.Line} = true diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 70243e8a..1b34405e 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -12,8 +12,11 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussLegendre; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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) @@ -40,8 +43,11 @@ function integral( f::F, plane::Meshes.Plane, rule::GaussKronrod; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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) @@ -58,8 +64,11 @@ function integral( f::F, plane::Meshes.Plane, rule::HAdaptiveCubature; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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) @@ -87,3 +96,9 @@ function integral( # Reapply units 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 fce24530..9008b4a8 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -12,8 +12,11 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussLegendre; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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) @@ -38,8 +41,11 @@ function integral( f::F, ray::Meshes.Ray, rule::GaussKronrod; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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)) @@ -53,8 +59,11 @@ function integral( f::F, ray::Meshes.Ray, rule::HAdaptiveCubature; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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)) @@ -78,3 +87,9 @@ function integral( # Reapply units return value .* integrandunits end + +################################################################################ +# jacobian +################################################################################ + +_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true 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, diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index e4aac959..9fb2af9a 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 = Analytical(), 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,11 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::GaussKronrod; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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] @@ -38,7 +42,14 @@ function integral( f::F, tetrahedron::Meshes.Tetrahedron, rule::HAdaptiveCubature; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} _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 d0b418bb..c63e715b 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -9,20 +9,15 @@ # 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.Ngon{3}, + triangle::Meshes.Triangle, rule::GaussLegendre; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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) @@ -53,18 +48,15 @@ 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.Ngon{3}, + triangle::Meshes.Triangle, rule::GaussKronrod; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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] @@ -74,19 +66,15 @@ 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.Ngon{3}, + triangle::Meshes.Triangle, rule::HAdaptiveCubature; + diff_method::DM = Analytical(), FP::Type{T} = Float64 -) where {F <: Function, T <: AbstractFloat} +) 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φ ) @@ -104,3 +92,9 @@ function integral( # Apply a linear domain-correction factor 0.5 ↦ area(triangle) return 2 * Meshes.area(triangle) .* ∫ end + +################################################################################ +# jacobian +################################################################################ + +_has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true diff --git a/src/utils.jl b/src/utils.jl index 14de6e6b..551ac1e1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,6 +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::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) + +# 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 ################################################################################ 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..921c8d13 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,4 +8,34 @@ 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 + 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 45453a50..dd1713e8 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -9,3 +9,38 @@ p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") @test MeshIntegrals._units(p) == u"m" 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 + + # AutoEnzyme + @test AutoEnzyme() isa AutoEnzyme +end