Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

v0.14-DEV: Consolidate nomenclature to IntegrationRule #82

Merged
merged 15 commits into from
Sep 24, 2024
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "MeshIntegrals"
uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
version = "0.13.5"
version = "0.14.0-DEV"

[deps]
CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
Expand Down
6 changes: 3 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,9 @@ 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, algorithm, FP=Float64)
integral(f, geometry, rule, FP=Float64)
```
Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration algorithm, e.g. `GaussKronrod()`.
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.

Expand All @@ -37,7 +37,7 @@ 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 `algorithm` can also be specified as a third-argument.
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)
Expand Down
2 changes: 1 addition & 1 deletion docs/src/api.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ MeshIntegrals.surfaceintegral
MeshIntegrals.volumeintegral
```

## Integration Algorithms
## Integration Rules

```@docs
MeshIntegrals.GaussKronrod
Expand Down
4 changes: 2 additions & 2 deletions docs/src/supportmatrix.md
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
# Support Matrix

While this library aims to support all possible integration algorithms and **Meshes.jl**
While this library aims to support all possible integration rules and **Meshes.jl**
geometry types, some combinations are ill-suited and some others are simplu not yet
implemented. The following Support Matrix aims to capture the current development state of
all geometry/algorithm combinations. Entries with a green check mark are fully supported
all geometry/rule combinations. Entries with a green check mark are fully supported
and have passing unit tests that provide some confidence they produce accurate results.

In general, Gauss-Kronrod integration rules are recommended (and the default) for geometries
Expand Down
8 changes: 5 additions & 3 deletions src/MeshIntegrals.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,14 @@ module MeshIntegrals
include("utils.jl")
export jacobian, derivative, unitdirection

include("integration_algorithms.jl")
export GaussKronrod, GaussLegendre, HAdaptiveCubature
include("integration_rules.jl")
export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature

include("integral.jl")
export integral

include("integral_aliases.jl")
export integral, lineintegral, surfaceintegral, volumeintegral
export lineintegral, surfaceintegral, volumeintegral

# Integration methods specialized for particular geometries
include("specializations/BezierCurve.jl")
Expand Down
62 changes: 31 additions & 31 deletions src/integral.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,17 +4,17 @@

"""
integral(f, geometry)
integral(f, geometry, algorithm)
integral(f, geometry, algorithm, FP)
integral(f, geometry, rule)
integral(f, geometry, rule, FP)

Numerically integrate a given function `f(::Point)` over the domain defined by
a `geometry` using a particular `integration algorithm` with floating point
a `geometry` using a particular numerical integration `rule` with floating point
precision of type `FP`.

# Arguments
- `f`: an integrand function with a method `f(::Meshes.Point)`
- `geometry`: some `Meshes.Geometry` that defines the integration domain
- `algorithm`: optionally, the `IntegrationAlgorithm` used for the 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)

Note that reducing `FP` below `Float64` will incur some loss of precision. By
Expand All @@ -23,7 +23,7 @@ contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision
"""
function integral end

# If only f and geometry are specified, select default algorithm
# If only f and geometry are specified, select default rule
function integral(
f::F,
geometry::G
Expand All @@ -33,14 +33,14 @@ function integral(
_integral(f, geometry, rule)
end

# with algorithm and T specified
# with rule and T specified
function integral(
f::F,
geometry::G,
settings::I,
rule::I,
FP::Type{T} = Float64
) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat}
_integral(f, geometry, settings, FP)
) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule, T<:AbstractFloat}
_integral(f, geometry, rule, FP)
end


Expand All @@ -52,31 +52,31 @@ end
function _integral(
f,
geometry,
settings::GaussKronrod,
rule::GaussKronrod,
FP::Type{T} = Float64
) where {T<:AbstractFloat}
# Run the appropriate integral type
Dim = Meshes.paramdim(geometry)
if Dim == 1
return _integral_1d(f, geometry, settings, FP)
elseif Dim == 2
return _integral_2d(f, geometry, settings, FP)
elseif Dim == 3
return _integral_3d(f, geometry, settings, FP)
N = Meshes.paramdim(geometry)
if N == 1
return _integral_1d(f, geometry, rule, FP)
elseif N == 2
return _integral_2d(f, geometry, rule, FP)
elseif N == 3
return _integral_3d(f, geometry, rule, FP)
end
end

# GaussLegendre
function _integral(
f,
geometry,
settings::GaussLegendre,
rule::GaussLegendre,
FP::Type{T} = Float64
) where {T<:AbstractFloat}
N = Meshes.paramdim(geometry)

# Get Gauss-Legendre nodes and weights for a region [-1,1]^N
xs, ws = _gausslegendre(FP, settings.n)
xs, ws = _gausslegendre(FP, rule.n)
weights = Iterators.product(ntuple(Returns(ws), N)...)
nodes = Iterators.product(ntuple(Returns(xs), N)...)

Expand All @@ -95,21 +95,21 @@ end
function _integral(
f,
geometry,
settings::HAdaptiveCubature,
rule::HAdaptiveCubature,
FP::Type{T} = Float64
) where {T<:AbstractFloat}
Dim = Meshes.paramdim(geometry)
N = Meshes.paramdim(geometry)

integrand(t) = f(geometry(t...)) * differential(geometry, t)

# HCubature doesn't support functions that output Unitful Quantity types
# Establish the units that are output by f
testpoint_parametriccoord = fill(FP(0.5),Dim)
testpoint_parametriccoord = fill(FP(0.5), N)
integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
# Create a wrapper that returns only the value component in those units
uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
# Integrate only the unitless values
value = HCubature.hcubature(uintegrand, zeros(FP,Dim), ones(FP,Dim); settings.kwargs...)[1]
value = HCubature.hcubature(uintegrand, zeros(FP,N), ones(FP,N); rule.kwargs...)[1]

# Reapply units
return value .* integrandunits
Expand All @@ -122,29 +122,29 @@ end
function _integral_1d(
f,
geometry,
settings::GaussKronrod,
rule::GaussKronrod,
FP::Type{T} = Float64
) where {T<:AbstractFloat}
integrand(t) = f(geometry(t)) * differential(geometry, [t])
return QuadGK.quadgk(integrand, FP(0), FP(1); settings.kwargs...)[1]
integrand(t) = f(geometry(t)) * differential(geometry, (t))
return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
end

function _integral_2d(
f,
geometry2d,
settings::GaussKronrod,
rule::GaussKronrod,
FP::Type{T} = Float64
) where {T<:AbstractFloat}
integrand(u,v) = f(geometry2d(u,v)) * differential(geometry2d, [u,v])
∫₁(v) = QuadGK.quadgk(u -> integrand(u,v), FP(0), FP(1); settings.kwargs...)[1]
return QuadGK.quadgk(v -> ∫₁(v), FP(0), FP(1); settings.kwargs...)[1]
integrand(u,v) = f(geometry2d(u,v)) * differential(geometry2d, (u,v))
∫₁(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

# Integrating volumes with GaussKronrod not supported by default
function _integral_3d(
f,
geometry,
settings::GaussKronrod,
rule::GaussKronrod,
FP::Type{T} = Float64
) where {T<:AbstractFloat}
error("Integrating this volume type with GaussKronrod not supported.")
Expand Down
Loading