Docstrings

allocate_buffer(
    fct,
    domain::AbstractDomain{D,T};
    embedded_cubature::EmbeddedCubature{D,T}=default_embedded_cubature(domain),
    norm=LinearAlgebra.norm
) where {D,T}

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.

integrate(
    fct,
    domain::AbstractDomain{D,T};
    embedded_cubature::EmbeddedCubature{D,T}=default_embedded_cubature(domain),
    subdiv_algo=default_subdivision(domain),
    buffer=nothing,
    norm=LinearAlgebra.norm,
    atol=zero(T),
    rtol=(atol > zero(T)) ? zero(T) : sqrt(eps(T)),
    maxsubdiv=2^(13 + D),
    callback=(I, E, nb_subdiv, buffer) -> nothing,
) where {D,T}

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 K. The return type K must support addition and multiplication by scalars of type 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

• embedded_cubature::EmbeddedCubature{D,T}=default_embedded_cubature(domain): the embedded cubature. Each supported domain has a default_embedded_cubature.
• subdiv_algo=default_subdivision(domain): the subdivision algorithm, each domain has a default_subdivision.
• buffer=nothing: optional heap used by the adaptive algorithm. Reusing a buffer from allocate_buffer can reduce allocations when calling integrate repeatedly.
• norm=LinearAlgebra.norm: norm used to estimate the error.
• atol=zero(T): absolute tolerance.
• rtol=(atol > zero(T)) ? zero(T) : sqrt(eps(T)): relative tolerance.
• 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, E ≤ rtol * norm(I), or nb_subdiv ≥ maxsubdiv.

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

default_embedded_cubature(domain::DOM) where {DOM<:AbstractDomain}

Return a default embedded cubature for the domains:

• dimension 1:
  • Segment: SEGMENT_GK15
• dimension 2:
  • Triangle: RadonLaurie()
  • Rectangle: SQUARE_CH25
• dimension 3:
  • Tetrahedron: GrundmannMoeller{3}(7, 5)
  • Cuboid: CUBE_BE65
• dimension D:
  • Simplex: GrundmannMoeller{D}(7, 5)
  • Orthotope: GenzMalik{D}()

default_subdivision(domain::DOM) where {DOM<:AbstractDomain}

Return the default algorithm to subdivide domain.

• dimension 1:
  • Segment: subdivide_segment
• dimension 2:
  • Triangle: subdivide_triangle
  • Rectangle: subdivide_rectangle
• dimension 3:
  • Tetrahedron: subdivide_tetrahedron
  • Cuboid: subdivide_cuboid
• dimension D:
  • Simplex: subdivide_simplex
  • Orthotope: subdivide_orthotope

abstract type AbstractDomain{D,T}

Abstract type for integration domains in D dimensions with element type T.

Type Parameters:

• D: dimension of the domain.
• T: element type of the domain.

Mandatory methods:

• map_from_reference

Useful (but non-mandatory) methods:

• reference_domain
• map_to_reference

const Cuboid{T} = Orthotope{3,T}

Alias for a 3-dimensional Orthotope of value type T.

Constructors:

• Cuboid(low_corner, high_corner)
• Cuboid{T}(low_corner, high_corner)

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}})

const Rectangle{T} = Orthotope{2,T}

Alias for a 2-dimensional Orthotope of element type T.

Constructors:

• Rectangle(low_corner, high_corner)
• Rectangle{T}(low_corner, high_corner)

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)

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}})

const Tetrahedron{T} = Simplex{3,T,4}

Alias for a 3-dimensional Simplex with 4 vertices of element type T.

Constructors:

• Tetrahedron(a, b, c, d)
• Tetrahedron{T}(a, b, c, d)

const Triangle{T} = Simplex{2,T,3}

Alias for a 2-dimensional Simplex with 3 vertices of value type T.

Constructors:

• Triangle(a, b, c)
• Triangle{T}(a, b, c)

dimension(::Type{<:AbstractDomain{D}}) where {D}

Return the dimension D encoded in the domain type.

Φ, μ = map_from_reference(domain::DOM) where {DOM<:AbstractDomain}

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

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.

reference_domain(::Type{<:AbstractDomain})

Return the reference domain for the given domain type. See reference_segment, reference_simplex, and reference_orthotope for more details.

reference_orthotope(::Val{D}, (::Type{T})=float(Int)) where {D,T}

Return the reference D-dimensional orthotope [0, 1]ᴰ with element type T.

reference_segment(T=float(Int))

Return the reference 1-dimensional segment [0, 1] with element type T.

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).

subdivide_cuboid(c::Cuboid)

Divide the cuboid c into 8 smaller cuboids by splitting each axis at its midpoint.

subdivide_orthotope(h::Orthotope)

Subdivide h into 2^D smaller orthotopes by bisecting each axis at its midpoint.

subdivide_rectangle(r::Rectangle)

Divide the rectangle r into 4 smaller rectangles by splitting each axis at its midpoint.

subdivide_segment(s::Segment)

Subdivide s into two smaller segments by splitting at its midpoint.

subdivide_simplex(s::Simplex{D,T,N}) where {D,T,N}

Subdivide a D-simplex into 2ᴰ simplices by using the SimpleS algorithm.

Implements the SimpleS algorithm in Algorithm 860: SimpleS – an extension of Freudenthal's simplex subdivision.

subdivide_simplex_scheme(::Val{D}) where {D}

Return the color schemes for subdividing a D-simplex into 2ᴰ simplices by using the SimpleS algorithm.

subdivide_tetrahedron(t::Tetrahedron)

Divide the tetrahedron t into 8 tetrahedra by connecting the midpoints of the edges.

subdivide_triangle(t::Triangle)

Divide the triangle t into 4 triangles by connecting the midpoints of the edges.

Berntsen-Espelid tabulated embedded cubature for Cuboid.

This rule uses 115 nodes with high order 11 and low order 9.

Berntsen-Espelid tabulated embedded cubature for Cuboid.

This rule uses 65 nodes with high order 9 and low order 7.

Gauss-Kronrod tabulated embedded cubature for Segment.

This rule uses 15 nodes with high order 23 and low order 13.

Gauss-Kronrod tabulated embedded cubature for Segment.

This rule uses 31 nodes with high order 47 and low order 29.

Gauss-Kronrod tabulated embedded cubature for Segment.

This rule uses 7 nodes with high order 11 and low order 5.

Cools-Haegemans tabulated embedded cubature for Rectangle.

This rule uses 21 nodes with high order 7 and low order 5.

Cools-Haegemans tabulated embedded cubature for Rectangle.

This rule uses 25 nodes with high order 9 and low order 7.

abstract type AbstractRule{DOM<:AbstractDomain}

Abstract type for a cubature rule on a domain DOM.

Type Parameters:

• DOM: reference_domain(DOM) gives the reference domain on which the embedded cubature is assumed to be defined.

Mandatory methods:

• embedded_cubature

Useful (but non-mandatory) methods:

• orders

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)

(ec::EmbeddedCubature{D,T})(
    fct, domain::AbstractDomain{D,T}, norm=LinearAlgebra.norm
) where {D,T}

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.

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.

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.

struct RadonLaurie <: AbstractRule{Simplex{2}}

Embedded cubature rule for a 2-simplex (a.k.a a triangle) of high order 8 and low order 5.

@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

embedded_cubature(ar::AbstractRule, T=float(Int))

embedded_cubature(nodes, weights_high, weights_low, T=float(Int))

Construct an embedded cubature with element type T.

The constructor can be called from a subtype of AbstractRule, or from explicit nodes, weights_high, and weights_low data. Available rule types include:

• TabulatedEmbeddedCubature
• RadonLaurie
• GrundmannMoeller
• GenzMalik

orders(rule::AR) where {AR<:AbstractRule}

Return (order_high, order_low) for rule.

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.