Docstrings

Integration API

HAdaptiveIntegration.allocate_buffer — Method
allocate_buffer(
    fct::FCT,
    domain::DOM;
    rule::EmbeddedCubature{D}=default_rule(domain),
    norm=LinearAlgebra.norm,
) where {D,FCT,DOM<:AbstractDomain{D}}

Allocate and return a heap buffer compatible with integrate.

Passing this buffer through the buffer keyword can reduce memory allocations when integrate is called repeatedly with compatible domain and value types.

source
HAdaptiveIntegration.integrate — Method
integrate(
    fct,
    domain::AbstractDomain{D};
    rule::EmbeddedCubature{D}=default_rule(domain),
    subdiv_algo=default_subdivision(domain),
    norm=LinearAlgebra.norm,
    buffer=nothing,
    atol=nothing,
    rtol=nothing,
    maxsubdiv=2^(13 + D),
    callback=(I, E, nb_subdiv, buffer) -> nothing,
) where {D}

Adaptively integrate fct over domain.

Return (I, E) where I is the integral estimate and E is an error estimate computed from an embedded cubature pair.

Arguments

  • fct: a function that maps SVector{D,T} to a value in K. The return type K must support addition and multiplication by scalars of type typeof(one(T)).
  • domain::AbstractDomain{D,T}: the integration domain. Currently, we support Segment, Triangle, Rectangle, Tetrahedron, Cuboid, D-dimensional Simplex, and D-dimensional Orthotope.

Optional arguments

  • rule::EmbeddedCubature{D,T}=default_rule(domain): the embedded cubature rule. Each supported domain has a default_rule.
  • subdiv_algo=default_subdivision(domain): the subdivision algorithm, each domain has a default_subdivision.
  • norm=LinearAlgebra.norm: norm used to estimate the error.
  • buffer=nothing: optional heap used by the adaptive algorithm. Reusing a buffer from allocate_buffer can reduce allocations when calling integrate repeatedly.
  • atol=nothing: absolute tolerance, if nothing is passed, it will be set to zero(E).
  • rtol=nothing: relative tolerance, if nothing is passed, it will be set to iszero(atol) ? sqrt(eps(one(E))) : zero(one(E)).
  • maxsubdiv=2^(13 + D): maximum number of subdivisions.
  • callback=(I, E, nb_subdiv, buffer) -> nothing: a callback function called for each estimate of I and E, including the initial estimate (nb_subdiv=0) and after each subdivision. The callback receives the current integral I, error estimate E, number of subdivisions nb_subdiv, and buffer.

Notes

  • Iteration stops when E ≤ atol or E ≤ rtol * norm(I) or nb_subdiv == maxsubdiv.
  • For this function to be thread safe, you need three conditions:
    1. fct must be thread safe.
    2. If a buffer is provided via the use of allocate_buffer, you need to create a new buffer for each thread. The default is thread safe.
    3. If a callback function is provided, it must be thread safe. The default is thread safe.
source
HAdaptiveIntegration.resum — Method
resum(buffer; norm=LinearAlgebra.norm)

Re-sum integral and error contributions stored in buffer.

This is more expensive than a plain sum, but it uses the Kahan-Babuška-Neumaier [1] summation algorithm to reduce floating-point round-off error.

[1] Klein, A. A Generalized Kahan-Babuška-Summation-Algorithm. Computing 76, 279-293 (2006). https://doi.org/10.1007/s00607-005-0139-x

source

Domains

HAdaptiveIntegration.Domain.Orthotope — Type
struct Orthotope{D,T} <: AbstractDomain{D,T}

Axis-aligned orthotope in D dimensions, with element type T, defined by two points low_corner and high_corner.

Fields:

  • corners::SVector{2,SVector{D,T}}: corners[1] is the low corner and corners[2] is the high corner.

Invariants (checked by outer constructors):

  • corners[1] .≤ corners[2]

Constructors:

  • Orthotope(low_corner, high_corner)
  • Orthotope{T}(low_corner, high_corner)
  • Orthotope(corners::SVector{2,SVector{D,T}})
source
HAdaptiveIntegration.Domain.Segment — Type
struct Segment{T} <: AbstractDomain{1,T}

One-dimensional segment with element type T, defined by xmin and xmax.

Fields:

  • xmin::T: lower endpoint.
  • xmax::T: upper endpoint.

Invariants (checked by outer constructors):

  • xmin < xmax

Constructors:

  • Segment(xmin, xmax)
  • Segment{T}(xmin, xmax)
source
HAdaptiveIntegration.Domain.Simplex — Type
struct Simplex{D,T,N} <: AbstractDomain{D,T}

A simplex in D dimensions with N=D+1 vertices of element type T.

Fields:

  • vertices::SVector{N,SVector{D,T}}: vertices of the simplex.

Invariants (not check at construction):

  • N = D+1

Constructors:

  • Simplex(vertices...)
  • Simplex{T}(vertices...)
  • Simplex(vertices::SVector{N,SVector{D,T}})
source
HAdaptiveIntegration.Domain.map_from_reference — Function
Φ, μ = map_from_reference(domain::DOM) where {DOM<:AbstractDomain}

Return (Φ, μ), where Φ is an anonymous function that maps the reference domain to the physical domain domain, and μ is the absolute value of the Jacobian determinant of Φ.

source
HAdaptiveIntegration.Domain.map_to_reference — Function
Ψ = map_to_reference(domain::DOM) where {DOM<:AbstractDomain}

Return Ψ an anonymous function that maps the physical domain to the reference domain.

Constraints:

  • For Segment, must have xmax ≠ xmin.
  • For Simplex{D}, the vertices must form a valid D-dimensional simplex with non-zero volume.
  • For Orthotope, must have high_corner ≠ low_corner.
source
HAdaptiveIntegration.Domain.reference_simplex — Method
reference_simplex(::Val{D}, (::Type{T})=float(Int)) where {D,T}

Return the reference D-dimensional simplex with element type T, which is the convex hull of the N=D+1 points (0,...,0), (1,0,...,0), (0,1,0,...,0), ..., (0,...,0,1).

source

Rules

HAdaptiveIntegration.Rule.EmbeddedCubature — Type
struct EmbeddedCubature{D,T}

An embedded cubature rule consisting of a high order cubature rule with H nodes and a low order cubature rule with L nodes. Note that the low order cubature uses nodes[1:L] as its nodes. The cubature nodes and weights are assumed to be defined on the reference domain.

Fields:

  • nodes::Vector{SVector{D,T}}: the cubature nodes.
  • weights_high::Vector{T}: the cubature weights for the high order cubature.
  • weights_low::Vector{T}: the cubature weights for the low order cubature.

Invariants (check at construction):

  • length(nodes) == length(weights_high)
  • length(weights_high) ≥ length(weights_low)
source
HAdaptiveIntegration.Rule.EmbeddedCubature — Method
(ec::EmbeddedCubature{D})(
    fct, domain::AbstractDomain{D}, norm=LinearAlgebra.norm
) where {D}

Return I_high and norm(I_high - I_low) where I_high and I_low are the result of the high order cubature and the low order cubature on domain. The function fct must take a SVector{D,T} to a return type F, and F must support scalar multiplication and addition. Note that there is no check, beyond compatibility of dimension and type, that the embedded cubature matches the intended domain geometry.

source
HAdaptiveIntegration.Rule.GenzMalik — Type
struct GenzMalik{D} <: AbstractRule{Orthotope{D}}

Embedded cubature rule for a D-orthotope of high order 7 and low order 5.

Type Parameters:

  • D: The dimension of the orthotope.
source
HAdaptiveIntegration.Rule.GrundmannMoeller — Type
struct GrundmannMoeller{D} <: AbstractRule{Simplex{D}}

Embedded cubature rule for a D-simplex.

Type Parameters:

  • D: The dimension of the simplex.

Fields:

  • order_high::Int: the high order of the rule.
  • order_low::Int: the low order of the rule.

Invariants (check at construction):

  • order_high and order_low must be odd.
  • must have order_high > order_low ≥ 1.
source
HAdaptiveIntegration.Rule.TabulatedEmbeddedCubature — Type
@kwdef struct TabulatedEmbeddedCubature{DOM<:AbstractDomain} <: AbstractRule{DOM}

An embedded cubature rule consisting of a high order cubature rule and a low order cubature rule. The low-order cubature uses nodes[1:L], where L = length(weights_low). Nodes and weights are defined on the reference domain (use reference_domain to inspect it).

Fields:

  • description::String: description of the embedded cubature.
  • reference::String: source of the tabulated values.
  • order_high::Int: order of the high order cubature.
  • order_low::Int: order of the low order cubature.
  • precision::Int: number of significant digits on the node and weight values, 10^-precision gives the approximate relative precision.
  • nodes::Vector{Vector{String}}: the cubature nodes.
  • weights_high::Vector{String}: the cubature weights for the high order cubature.
  • weights_low::Vector{String}: the cubature weights for the low order cubature.

Invariants (check at construction):

  • length(nodes) == length(weights_high)
  • length(weights_high) ≥ length(weights_low)
  • order_high ≥ order_low ≥ 0
  • precision ≥ 0
source

Extensions

HAdaptiveIntegration.increase_precision — Function
increase_precision(
    tec::TabulatedEmbeddedCubature,
    ::Type{T};
    x_atol=10 * eps(T),
    f_atol=10 * eps(T),
    maxiter::Int=16,
)

Increase the precision of a TabulatedEmbeddedCubature to match the precision of the floating-point type T. This method is provided by an extension, so ForwardDiff must be loaded.

Arguments

  • tec::TabulatedEmbeddedCubature: the tabulated embedded cubature to increase the precision of.
  • ::Type{T}: target floating-point type.

Optional arguments

  • x_atol=10 * eps(T): the absolute tolerance for the change in the variables (nodes and weights) between Newton iterations.
  • f_atol=10 * eps(T): the absolute tolerance for the change in the function values (integral constraints) between Newton iterations.
  • maxiter::Int=16: the maximum number of Newton iterations to perform.
source