References

const ACOUSTIC_SIZE_THRESHOLD

Boxes with a size larger than ACOUSTIC_SIZE_THRESHOLD times the wavelength are considered "high-frequency", and their cone-domain size will be refined.

const SAME_POINT_TOLERANCE = 1e-12

For singular kernels, two points are considered the same if their distance is less than SAME_POINT_TOLERANCE.

abstract type AbstractPDE{N}

A partial differential equation in dimension N. AbstractPDE types are used to define PDEKernels.

abstract type AbstractSplitter

An AbstractSplitter is used to split a ClusterTree. The interface requires the following methods:

• should_split(clt,splitter) : return a Bool determining if the ClusterTree should be further divided
• split!(clt,splitter) : perform the splitting of the ClusterTree handling the necessary data sorting.

See DyadicSplitter for an example of an implementation.

struct AdmissibilityCondition

Functor used to determine if the interaction between a target and a source can be approximated using an interpolation scheme.

A target and source are admissible under this admissiblility condition if the target box lies farther than r*η away, where r is the radius of the source box and η >= 1 is an adjustable parameter. In dimension N, the default choice is η = √N.

mutable struct ClusterTree{T,S,D}

Tree structure used to cluster elements of type V = eltype(T) into containers of type S. The method center(::V)::SVector is required for the clustering algorithms. An additional data field of type D can be associated with each node to store node-specific information (it defaults to D=Nothing).

Fields:

• _elements::T : vector containing the sorted elements.
• container::S : container for the elements in the current node.
• index_range::UnitRange{Int} : indices of elements contained in the current node.
• loc2glob::Vector{Int} : permutation from the local indexing system to the original (global) indexing system used as input in the construction of the tree.
• children::Vector{ClusterTree{N,T,D}}
• parent::ClusterTree{N,T,D}
• data::D : generic data field of type D.

ClusterTree(elements,splitter;[copy_elements=true])
ClusterTree{D}(points,splitter;[copy_elements=true])

Construct a ClusterTree from the given elements using the splitting strategy encoded in splitter. If copy_elements is set to false, the elements argument are directly stored in the ClusterTree and are permuted during the tree construction.

struct ConeDomainSize

Use ConeDomainSize(k,ds) to create a function that computes the appropriate meshsize for a given ClusterTree based on the wavenumber k and the reference meshsize ds.

Calling f(node) will return the meshsize for the given node, where f is a ConeDomainSize object.

struct DyadicMaxDepthSplitter <: AbstractSplitter

Similar to DyadicSplitter, but splits nodes until a maximum depth is reached.

See also: AbstractSplitter

struct DyadicMinimalSplitter <: AbstractSplitter

Similar to DyadicSplitter, but the boundin boxes are shrank to the minimal axis-aligned boxes at the end.

See also: AbstractSplitter

struct DyadicSplitter <: AbstractSplitter

Used to split an N dimensional ClusterTree into 2^N children until at most nmax points are contained in node.

See also: AbstractSplitter

struct Helmholtz{N, T}

Helmholtz equation in N dimension: Δu + k²u= 0. k is the wavenumber of numeric type T.

struct HyperRectangle{N,T}

Axis-aligned hyperrectangle in N dimensions given by low_corner::SVector{N,T} and high_corner::SVector{N,T}. By default, T is Float64.

HyperRectangle(lc,hc)

Construct an axis-aligned hyperrectangle with low-corner lc and high-corner hc.

It is recommended that lc and hc be of type SVector or Tuple to avoid type instabilities.

julia> IFGF.HyperRectangle((0.0, 0.0), (1.0, 2.0))
IFGF.HyperRectangle{2, Float64}([0.0, 0.0], [1.0, 2.0])

struct IFGFOp{T} <: AbstractMatrix{T}

High-level structure representing a linear operator with elements of type T and i,j entry given by K(X[i],Y[j]), where K is the underlying kernel function, and X and Y are the elements of the target_tree and source_tree, respectively.

Fields:

• kernel : underlying kernel function
• target_tree : tree structure used for clustering the target points
• source_tree : tree structure used for clustering the source points
• p::NTuple: number of chebyshev poins per dimension used in the interpolation of far interactions
• h::NTuple: meshsize in parameter (per dimension) of the largest cones

struct Laplace{N}

Laplace equation in N dimension: Δu = 0.

struct MulBuffer{T}

Low-level structure used to hold the interpolation coefficients and values of the current and next level cone interpolants.

The type paramter T is the type of the interpolant.

abstract type PDEKernel{Op}

A kernel function associated a Green function (or derivatives of) a partial differential equation.

const SourceTree{T,S,D}

Type alias for a ClusterTree with a data field D of type SourceTreeData.

mutable struct SourceTreeData{N,T,S}

Struct used to store data associated with a node of the SourceTree. The parametric type N represents the ambient dimension, T represents the primitive type used to represent the points (e.g. Float32 or Float64 for single/double precision), and S is the type of node in the near and far field.

Fields:

• msh: uniform cartesian mesh of the interpolation space, indexable by a CartesianIndex
• conedatadict: dictionary mapping a cone tag I::CartesianIndex to the index range of the interpolation points owned by that cone
• farlist: nodes on the target tree that can be evaluated through interpolation
• nearlist: nodes on the target tree that must be evaluated by a direct sum.

struct Stokes{N}

Stokes equation in N dimension.

const TargetTree{T,S,D}

Type alias for a ClusterTree with a data field D of type TargetTreeData.

struct TargetTreeData

Struct used to store data associated with a node of the TargetTree. Currently there is nothing there, but maybe one day there will be a downward pass in the target tree, and it will need some fields to store data.

struct UniformCartesianMesh{N,T} <: AbstractMesh{N,T}

An N-dimensional cartesian grid given as the tensor-product of N one-dimensional LinRange{T} grids.

Iterating over a UniformCartesianMesh generates the elements which compose the mesh; i.e. the HyperRectangle cells.

UniformCartesianMesh(domain::HyperRectangle,sz::NTuple)
UniformCartesianMesh(domain::HyperRectangle;step::NTuple)

Construct a uniform UniformCartesianMesh with sz[d] elements along dimension d. If the kwarg step is passed, construct a UniformCartesianMesh with elements of approximate size step.

split(rec::HyperRectangle[, axis::Int, place])

Split a hyperrectangle in two along the axis direction at the position place. Returns a tuple with the two resulting hyperrectangles.

When no place is given, defaults to splitting in the middle of the axis.

When no axis and no place are given, defaults to splitting along the largest axis.

union(r1::HyperRectangle,r2::HyperRectangle)

Smallest HyperRectangle that contains both r1 and r2.

_addnewcone!(node::SourceTree,I::CartesianIndex[, lock])

If node does not contain cone I, add it to the conedict. Does not set the indices of the interpolation data owned by the cone.

If lock is passed, the operation is thread-safe.

_binary_split!(cluster::ClusterTree,dir,pos;parentcluster=cluster)
_binary_split!(cluster::ClusterTree,f;parentcluster=cluster)

Split a ClusterTree into two, sorting all elements in the process. For each resulting child assign child.parent=parentcluster.

Passing a dir and pos arguments splits the container of node along direction dir at position pos, then sorts all points into the resulting left/right nodes.

If passed a predicate f, each point is sorted according to whether f(x) returns true (point sorted on the left node) or false (point sorted on the right node). At the end a minimal HyperRectangle containing all left/right points is created.

_compute_cone_list!(partition, p, ds_func)

For each node in partition, compute the cones which are necessary to cover both the target points on its farlist as well as the interpolation points on its parent's cones. p denotes the number of interpolation points per dimension, and ds_func(node) gives the meshsize of the domain in interpolation-space for node.

_compute_interaction_list!(source, target, adm)

For each node in source, compute its nearlist and farlist. The nearlist contains all nodes in target for which the interaction must be computed using the direct summation. The farlist includes all nodes in target_tree for which an accurate interpolant in available in the source node.

The adm function is called through adm(target,source) to determine if target is in the farlist of source. When that is not the case we recurse on their children unless both are a leaves, in which case target is added to the nearlist of source.

_density_type_from_kernel_type(T)

Helper function to compute the expected density V type associated with a kernel type T. For T<:Number, simply return T. For T<:SMatrix, return an SVector with elements of the type eltype(T) and of length size(T,2) so that multiplying an element of type T with an element of type V makes sense.

_forward_map!(out, L, c, v, grad)

Execute the plan of _plan_forward_map and store the result in out.

_init_msh!(node,domain,ds)

Create a uniform cartesian mesh of the interpolation domain of node with element spacing ds.

_plan_forward_map(pde, targets, sources, c, v, grad; kwargs...)

Assemble the IFGFOp for the forward map given by the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^N G(\boldsymbol{x}_i - \boldsymbol{y}_j) c_j + \gamma_{1,\boldsymbol{y}} G(\boldsymbol{x}_i - \boldsymbol{y}_j) v_j,\]

where $G$ is a Green's function of the PDE, $\gamma_{1,\boldsymbol{y}}$ its generalized Neumann trace respect to the $\boldsymbol{y}$ variable, $c_j$ denotes the value of x[j], and $v_j$ denotes the value of v[j]. If grad=true, the gradient of the expression is computed instead. Either c or v can be nothing, in which case the corresponding term is skipped in the sum.

This function handles the logic needed to determine the appropriate kernel function, as well values of p (number of interpolation points) and h (cone meshsize), and then forwards the keyword arguments to assemble_ifgf.

_unsafe_wrap_vector_of_sarray(A::Array)

Wrap the memory in A so that its first n-1 dimensions are interpreted as an SArray, returning a Vector{<:SArray} type.

This function uses unsafe_wrap, and therefore is unsafe for the same reasons unsafe_wrap is.

assemble_ifgf(K, X, Y; tol, nmax=250)
assemble_ifgf(K, X, Y; p, h, nmax=250)

Construct an approximation of the matrix [K(x,y) for x ∈ X, y ∈ Y] using the interpolated factored Green function method. The parameter tol specifies the desired (relative) error in Frobenius norm of the approximation, and nmax gives the maximum number of points in a leaf of the dyadic trees used.

The kernel K must implement wavenumber(::typeof(K)) --> Number to return the kernel's (typical) wavenumber; for non-oscillatory kernels this is zero.

When tol is passed, the number of interpolation points per dimension, p, and the meshsize h is estimated based on some heuristic. If you know in advance a good value of p and h for you function K and desired tolerance, it is suggested you pass that as a keyword argument instead.

bounding_box(els[, cube=false])

Compute a HyperRectangle that contains the center of all elements in els. If cube==true, the result is a square/cube instead.

Works for any type that implements center.

cart2interp(x,c,h)
cart2interp(x,source::SourceTree)

Map cartesian coordinates x into interpolation coordinates s centered at c and scaled by h. Interpolation coordinates are similar to hyperspherical coordinates (e.g. polar in two-dimensions and spherical in three-dimensions), except that the radial coordinates is replaced by s=h/r.

If passed a SourceTree, use its center and radius as c and h respectively.

center(Ω)

Center of the smallest ball containing Ω.

centered_factor(K,x,Y::SourceTree)

Return a representative scalar value of K(x,y) for y ∈ Y. This is used to help the interpolation of I(x) = ∑ᵢ K(x,yᵢ) where yᵢ ∈ Y and points x in the far field of Y. Defaults to K(x,yc) where yc is the center of Y.

cheb_error_estimate(coefs::AbstractArray{T,N},dim)

Given an N dimensional array of coefficients , estimate the relative error commited by the Chebyshev interpolant along dimension dim.

chebeval(coefs,x,rec[,sz])

Evaluate the Chebyshev polynomial defined on rec with coefficients given by coefs at the point x. If the size of coefs is known statically, its size can be passed as a Val using the sz argument.

chebnodes(n,[domain])
chebnodes(p::NTuple,[domain])
chebnodes(::Val{P},[domain])

Return the n Chebyshev nodes of the second kind on the interval [-1,1].

If passed a tuple p, return the tensor product of one-dimensinal nodes of size p.

If passed a Val{P}, statically compute the tensor product nodes.

Finally, passing a domain::HyperRectangle shifts the interpolation nodes to it.

lb = (0.,0.)
ub = (3.,2.)
domain = HyperRectangle(lb,ub)
chebnodes((10,5),domain) # generate `10×5` Chebyshev nodes on the `[0,3]×[0,2]` rectangle.

chebtransform_fftw!(vals::AbstractArray[, plan])

Given a function values on the chebnodes(p), where p=size(vals), compute the Chebyshev coefficients in-place. The underlying implementation uses the FFTW library. A pre-computed plan::FFTW.FFTWPlan should be passed for efficiency.

chebtransform_native!(vals)

Compute the Chebyshev coefficients given the values of a function on the Chebyshev nodes of the second kind (i.e. the extrema of the Chebyshev polynomials Tₙ). The result is written in-place.

cone_index(x::SVector,tree::SourceTree)

Given a point x in cartesian coordinates, return the index I of the cone in source to which x belongs, as well as its parametric coordinate s.

See also: cart2interp

container(clt::ClusterTree)

Return the object enclosing all the elements of the clt.

container_type(clt::ClusterTree)

Type used to enclose the elements of clt.

coords(x)

Return an SVector with the cartesian coordinates associated to a geometrical object x.

decrement_index(I::CartesianIndex,k[,n=1])

Equivalent to increment_index(I,k,-n)

depth(tree,acc=0)

Recursive function to compute the depth of node in a a tree-like structure.

Overload this function if your structure has a more efficient way to compute depth (e.g. if it stores it in a field).

diameter(Ω)

Largest distance between x and y for x,y ∈ Ω.

distance(X::ClusterTree, Y::ClusterTree)

Distance between the containers of X and Y.

distance(Ω₁,Ω₂)

Minimal Euclidean distance between a point x ∈ Ω₁ and y ∈ Ω₂.

distance(x::SVector,r::HyperRectangle)

The (minimal) Euclidean distance between the point x and any point y ∈ r.

element_index(s::SVector,m::UniformCartesianMesh)

Given a point s ∈ m, return the index I of the element in m containing s.

element_index_for_point(p::SVector,m::UniformCartesianMesh)

Given a point p, return the index I of the element in m containing p.

element_type(clt::ClusterTree)

Type of elements sorted in clt.

elements(clt::ClusterTree)

Iterable list of the elements inside clt.

float_type(x)

Floating type used to represent the underlying structure (e.g. Float32, Float64).

glob2loc(clt::ClusterTree)

The inverse of loc2glob(clt).

helmholtz3d(k, targets, sources; charges = nothing, dipvecs = nothing, grad =
false, kwargs...)

Compute the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^n \frac{c_j e^{i k || \boldsymbol{x}_i - \boldsymbol{y}_j ||}}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||} + \nabla_{\boldsymbol{y}}^\top \left( \frac{e^{i k || \boldsymbol{x}_i - \boldsymbol{y}_j ||}}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||} \right) v_j,\]

where $c_j$ are the charges and v_j are the dipvecs.

Input:

• k: the wavenumber.
• targets::Matrix{T}: 3 x n matrix of target points.
• sources::Matrix{T}: 3 x m matrix of source points.
• charges::Vector{Complex{T}}: vector of n charges.
• dipvecs::Matrix{Complex{T}}: 3 x m matrix of dipole vectors.
• grad::Bool: if true, the gradient is computed instead
• kwargs...: additional keyword arguments passed to assemble_ifgf.

Output

• u: the potential at the target points if grad=false, or the gradient as a 3 x n matrix if grad=true.

increment_index(I::CartesianIndex,k[,n=1])

Increment I by n along the dimension k. This is equivalent to I += n*eₖ, where eₖ is a vector with with 1 at the k-th coordinate and zeros elsewhere.

index_range(clt::ClusterTree)

Indices of elements in root_elements(clt) which lie inside clt.

interp2cart(s,xc,h)
interp2cart(s,source::SourceTree)

Map interpolation coordinate s into the corresponding cartesian coordinates.

Interpolation coordinates are similar to polar (in 2d) or spherical (in 3d) coordinates centered at the xc, except that the radial coordinate r is replaced by s = h/r.

If passed a SourceTree as second argument, call interp2cart(s,xc,h) with xc the center of the SourceTree and h its radius.

interpolation_domain(source::SourceTree)

Compute the domain (as a HyperRectangle) in interpolation space for which source must provide a covering through its cone interpolants.

inv_centered_factor(K,x,Y)

Inverse of centered_factor.

laplace3d!(out, L::IFGFOp; charges = nothing, dipvecs = nothing, grad =
false)

In-place version of laplace3d.

laplace3d(targets, sources; charges = nothing, dipvecs = nothing, grad =
false, kwargs...)

Compute the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^n \frac{c_j}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||} - \frac{(\boldsymbol{x}_i - \boldsymbol{y}_j)^\top v_j}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||^3} ,\]

where $c_j$ are the charges and v_j are the dipvecs.

Input:

• targets::Matrix{T}: 3 x n matrix of target points.
• sources::Matrix{T}: 3 x m matrix of source points.
• charges::Vector{T}: vector of n charges.
• dipvecs::Matrix{T}: 3 x m matrix of dipole vectors.
• grad::Bool: if true, the gradient is computed instead
• kwargs...: additional keyword arguments passed to assemble_ifgf (e.g. tol)

Output

• u: the potential at the target points if grad=false, or the gradient as a 3 x n matrix if grad=true.

loc2glob(clt::ClusterTree)

The permutation from the (local) indexing system of the elements of the clt to the (global) indexes used upon the construction of the tree.

near_interaction!(C,K,X,Y,σ,I,J)

Compute C[i] <-- C[i] + ∑ⱼ K(X[i],Y[j])*σ[j] for i ∈ I, j ∈ J. The default implementation simply does a double loop over all i ∈ I and j ∈ J: override this method for the type of your own custom kernel K if you have a fast way of computing the full product in-place (e.g. using SIMD instructions).

partition_by_depth(tree)

Given a tree, return a partition vector whose i-th entry stores all the nodes in tree with depth=i-1. Empty nodes are not added to the partition.

plan_laplace3d(targets, sources; charges = nothing, dipvecs = nothing, grad =
false, kwargs...)

Similar to laplace3d, but returns a plan in the form of an IFGFOpobject instead of the result. Callinglaplace3dwithplan` as the first argument will then compute the result.

radius(Ω)

Half the diameter.

real_and_imag(v)

Return a view over the real and imaginary parts of a complex vector v.

julia> v = [1.0 + 2.0im, 3.0 + 4.0im];

julia> vr,vi = IFGF.real_and_imag(v)
([1.0, 3.0], [2.0, 4.0])

return_type(f[,args...])

The type returned by f(args...), where args is a tuple of types. Falls back to Base.promote_op by default.

A functors of type T with a known return type should extend return_type(::T,args...) to avoid relying on promote_op.

root_elements(clt::ClusterTree)

The elements contained in the root of the tree to which clt belongs.

section(rec::HyperRectangle{D},ax)

Return the dimension D-1 HyperRectangle obtained by deleting the coordinates at the ax dimension.

share_interp_data(mode=true)

If set to true, neighboring cone interpolants will share their interpolation data (which are, in this case, computed only once).

should_split(clt::ClusterTree, depth, splitter::AbstractSplitter)

Determine whether or not a ClusterTree should be further divided.

split!(clt::ClusterTree,splitter::AbstractSplitter)

Divide clt using the strategy implemented by splitter. This function is reponsible of assigning the children and parent fields, as well as of permuting the data of clt.

stokes3d(targets, sources; stoklet = nothing, strslet = nothing, strsvec =
nothing, grad = false, kwargs...)

Compute the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^n \left( \frac{I}{d} + \frac{\boldsymbol{r}_{ij}\boldsymbol{r}_{ij}^\top}{2d_{ij}^3} \right) \boldsymbol{\sigma}_j + \frac{3(\boldsymbol{\nu}_j \cdot \boldsymbol{r}_{ij})(\boldsymbol{\mu}_j\cdot\boldsymbol{r}_{ij}) \boldsymbol{r}_{ij}}{d_{ij}^5},\]

where $\boldsymbol{r}_{ij} = |\boldsymbol{x}_i - \boldsymbol{y}_j|$, $d_{ij} = |\boldsymbol{r_{ij}}|$, $\boldsymbol{\sigma}_j$ are the Stokeslet strengths, $\boldsymbol{\mu}_j$ are the stresslet strengths, and $\boldsymbol{\nu}_j$ are the stresslet orientations.

Input:

• targets::Matrix{T}: 3 x n matrix of target points.
• sources::Matrix{T}: 3 x m matrix of source points.
• stoklet::Matrix{T}: 3 × n matrix of stokeslets
• strslet::Matrix{T}: 6 x m matrix of stresslet densities (indices 1:3) and orientations (indices 4:6)
• grad::Bool: if true, the gradient is computed instead
• kwargs...: additional keyword arguments passed to assemble_ifgf.

Output

• u: 3 x n matrix giving the velocity at the target points if grad=false, or 3 x 3 × n matrix representing the gradient of u if grad=true.

svector(f,n)

Create an SVector of length n, computing each element as f(i), where i ∈ 1:n.

transfer_factor(K,x,Y::SourceTree)

Ratio between the centered factor of Y and the centered factor of its parent.

use_fftw(mode=true)

Use the FFTW library for computing the Chebyshev transform.

use_minimal_conedomain(mode=true)

If set to true, the minimal domain for the cone interpolants will be used by calling interpolation_domain during the contruction of the active cones.

While setting this to true reduces the memory footprint by creating fewer cones, it also increases the time to assemble the IFGFOperator.

usethreads(mode=true)

Enable/disable thread usage globally for the IFGF module.

wavenumber(K::Function)

For oscillatory kernels, return the characteristic wavenumber (i.e 2π divided by he wavelength). For non-oscillatory kernels, return 0.

mul!(y,A::IFGFOp,x,a,b;global_index=true)

Setting global_index=false uses the internal indexing system of the IFGFOp, and therefore does not permute the vectors x and y to perform the multiplication.