Docstrings
HAdaptiveIntegration.Domain.AbstractDomain — Type
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:
Useful (but non-mandatory) methods:
HAdaptiveIntegration.Domain.Cuboid — Type
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)
HAdaptiveIntegration.Domain.Orthotope — Type
struct Orthotope{D,T} <: AbstractDomain{D,T}Axes-aligned Orthotope in D dimensions, with element type T, given by two points low_corner and high_corner. Note that, we must have low_corner .≤ high_corner.
Fields:
corners::SVector{2,SVector{D,T}}:corners[1]is the low corner andcorners[2]is the high corner.
Invariants (not check at construction):
corners[1] .≤ corners[2]
Constructors:
Orthotope(low_corner, high_corner)Orthotope{T}(low_corner, high_corner)Orthotope(corners::SVector{2,SVector{D,T}})
HAdaptiveIntegration.Domain.Rectangle — Type
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)
HAdaptiveIntegration.Domain.Segment — Type
struct Segment{T} <: AbstractDomain{1,T}Segment in 1 dimension, with element type T, given by xmin and xmax. Note that, we must have xmin ≤ xmax.
Fields:
xmin::T: the low corner.xmax::T: the high corner.
Invariants (not check at construction):
xmin ≤ xmax
Constructors:
Segment(xmin, xmax)
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}})
HAdaptiveIntegration.Domain.Tetrahedron — Type
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)
HAdaptiveIntegration.Domain.Triangle — Type
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)
HAdaptiveIntegration.Domain.combinations — Method
combinations(n::Int, k::Int)Helper function to generate all combinations of k elements from 1:n, similar to calling combinations(1:n, k) from Combinatorics.jl.
HAdaptiveIntegration.Domain.dimension — Method
dimension(::Type{<:AbstractDomain{D}}) where {D}Return the dimension D of the given domain DOM.
HAdaptiveIntegration.Domain.map_from_reference — Function
Φ, μ = map_from_reference(domain::DOM) where {DOM<:AbstractDomain}Return Φ and μ where Φ is an anonymous function that maps the reference domain to the physical domain domain and μ is the absolute value of the Jacobian's determinant of the map Φ.
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 havexmax > xmin. - For
Simplex{D}, the vertices must form a validD-dimensional simplex with non-zero volume). - For
Orthotope, must havehigh_corner .> low_corner.
HAdaptiveIntegration.Domain.reference_domain — Function
reference_domain(::Type{<:AbstractDomain})Return the reference domain for the given domain type. See reference_segment, reference_simplex, and reference_orthotope for more details.
HAdaptiveIntegration.Domain.reference_orthotope — Method
reference_orthotope(D::Int, T=float(Int))Return the reference D-dimensional orthotope [0, 1]ᴰ with element type T.
HAdaptiveIntegration.Domain.reference_segment — Method
reference_segment(T=float(Int))Return the reference 1-dimensional segment [0, 1] with element type T.
HAdaptiveIntegration.Domain.reference_simplex — Method
reference_simplex(D::Int, T=float(Int))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).
HAdaptiveIntegration.Domain.subdivide_cuboid — Method
subdivide_cuboid(c::Cuboid)Divide the cuboid c into 8 cuboid by connecting the center of the cuboid to the midpoints of the edges.
HAdaptiveIntegration.Domain.subdivide_orthotope — Method
subdivide_orthotope(h::Orthotope)Subdivide the D-orthotope h into 2ᴰ smaller orthotopes by splitting each dimension at its midpoint.
HAdaptiveIntegration.Domain.subdivide_rectangle — Method
subdivide_rectangle(r::Rectangle)Divide the rectangle r into 4 squares by connecting the center of the square to the midpoints of the edges.
HAdaptiveIntegration.Domain.subdivide_reference_orthotope — Method
subdivide_reference_orthotope(::Val{D}, ::Type{T}=float(Int))Like subdivide_orthotope, but operates on the reference orthotope. Since the output depends only on the dimension D, and the type T used to represent coordinates, this function is generated for each combination of D and T.
HAdaptiveIntegration.Domain.subdivide_reference_simplex — Method
subdivide_reference_simplex(::Val{D}, ::Type{T}=float(Int)) where {D,T}Like subdivide_simplex, but operates on the reference simplex. Since the output depends only on the dimension D, and the type T used to represent coordinates, this function is generated for each combination of D and T.
HAdaptiveIntegration.Domain.subdivide_segment — Method
subdivide_segment(s::Segment)Subdivide the 1-dimensional segment s into 2 smaller segments by splitting it at its midpoint.
HAdaptiveIntegration.Domain.subdivide_simplex — Method
subdivide_simplex(s::Simplex)Subdivide a D-simplex into 2ᴰ simplices by using the Freudenthal triangulation.
Implements the RedRefinementND algorithm in Simplicial grid refinement: on Freudenthal's algorithm and the optimal number of congruence classes.
HAdaptiveIntegration.Domain.subdivide_tetrahedron — Method
subdivide_tetrahedron(t::Tetrahedron)Divide the tetrahedron t into 8 tetrahedra by connecting the midpoints of the edges.
HAdaptiveIntegration.Domain.subdivide_triangle — Method
subdivide_triangle(t::Triangle)Divide the triangle t into 4 triangles by connecting the midpoints of the edges.
HAdaptiveIntegration.Rule.CUBE_BE115 — Constant
Berntsen-Espelid with 115 nodes.
HAdaptiveIntegration.Rule.CUBE_BE65 — Constant
Berntsen-Espelid with 65 nodes.
HAdaptiveIntegration.Rule.SEGMENT_GK15 — Constant
Gauss-Kronrod with 15 nodes.
HAdaptiveIntegration.Rule.SEGMENT_GK31 — Constant
Gauss-Kronrod with 31 nodes.
HAdaptiveIntegration.Rule.SEGMENT_GK7 — Constant
Gauss-Kronrod with 7 nodes.
HAdaptiveIntegration.Rule.SQUARE_CH21 — Constant
Cools-Haegemans with 21 nodes.
HAdaptiveIntegration.Rule.SQUARE_CH25 — Constant
Cools-Haegemans with 25 nodes.
HAdaptiveIntegration.Rule.AbstractRule — Type
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 assume to be set.
Mandatory methods:
Useful (but non-mandatory) methods:
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 for 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)
HAdaptiveIntegration.Rule.EmbeddedCubature — Method
(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 K, and K must support the multiplication by a scalar and the addition. Note that there is no check, beyond compatibility of dimension and type, that the embedded cubature is for the right domain.
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.
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_highandorder_lowmust be odd.- must have
order_high > order_low ≥ 1.
HAdaptiveIntegration.Rule.RadonLaurie — Type
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.
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. Note that the low order cubature uses nodes[1:L] as its nodes where L is the length of the weights_low. The cubature nodes and weights are assumed to be for the reference domain (use the reference_domain function to get the reference domain).
Fields:
description::String: description of the embedded cubature.reference::String: where the values are found.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^-precisiongive the relative precision of the values.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≥ 0precision ≥ 0
HAdaptiveIntegration.Rule.embedded_cubature — Method
embedded_cubature(ar::AbstractRule, T=float(Int))
embedded_cubature(nodes, weights_high, weights_low, T=float(Int))Construct the embedded cubature with element type T from a subtype of AbstractRule or from a vector of nodes and two vectors of weights for the high and low order cubature, with element type T. The list of AbstractRule's subtype are:
HAdaptiveIntegration.Rule.orders — Method
orders(rule::AR) where {AR<:AbstractRule}Return the high and low order of the embedded cubature rule.
HAdaptiveIntegration.allocate_buffer — Method
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 buffer that can be passed to the integrate function to improve performance by reducing memory allocations when integrate is called multiple times.
HAdaptiveIntegration.default_embedded_cubature — Method
default_embedded_cubature(domain::DOM) where {DOM<:AbstractDomain}Return a default embedded cubature for the domains:
- dimension 1:
- dimension 2:
- dimension 3:
Tetrahedron:GrundmannMoeller{3}(7, 5)Cuboid:CUBE_BE65
- dimension
D:Simplex:GrundmannMoeller{D}(7, 5)Orthotope:GenzMalik{D}()
HAdaptiveIntegration.default_subdivision — Method
default_subdivision(domain::DOM) where {DOM<:AbstractDomain}Return the default algorithm to subdivide domain.
- dimension 1:
- dimension 2:
- dimension 3:
- dimension
d:
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 is an extension, you need using ForwardDiff for using it.
Arguments
tec::TabulatedEmbeddedCubature: the tabulated embedded cubature to increase the precision of.::Type{T}: the floating point type to increase the precision to.
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.
HAdaptiveIntegration.integrate — Method
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}Return I and E where I is the integral of the function fct over domain and E is an error estimate.
Arguments
fct: a function that must take aSVector{D,T}to a return typeK, withKmust support the multiplication by a scalar of typeTand the addition.domain::AbstractDomain{D,T}: the integration domain. Currently, we supportSegment,Triangle,Rectangle,Tetrahedron,Cuboid,d-dimensionalSimplex, andd-dimensionalOrthotope.
Optional arguments
embedded_cubature::EmbeddedCubature{D,T}=default_embedded_cubature(domain): the embedded cubature, each supported domain has adefault_embedded_cubature.subdiv_algo=default_subdivision(domain): the subdivision algorithm, each domain has adefault_subdivision.buffer=nothing: heap use to do the adaptive algorithm, can be allocated usingallocate_buffer, which might result in performance gain if multiple call to integrate is perform.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 subdivision.callback=(I, E, nb_subdiv, buffer) -> nothing: a callback function called for each estimated value ofIandE, including the initial estimate (nb_subdiv=0) and after each subdivision. The callback receives the current integralI, error estimateE, number of subdivisionsnb_subdiv, andbuffer.
HAdaptiveIntegration.resum — Method
resum(buffer; norm=LinearAlgebra.norm)Re-sum the integral and error estimate from a provided buffer. This function is more expensive than a sum because it uses the Kahan-Babuška-Neumaier [1] summation algorithm to reduce numerical error due to floating-point numbers.
[1] Klein, A. A Generalized Kahan-Babuška-Summation-Algorithm. Computing 76, 279-293 (2006). https://doi.org/10.1007/s00607-005-0139-x