Docstrings

module Inti

Library for solving integral equations using Nyström methods.

const COMPRESSION_METHODS = [:none, :hmatrix, :fmm]

Available compression methods for the dense linear operators in Inti.

const CORRECTION_METHODS = [:none, :dim, :adaptive]

Available correction methods for the singular and nearly-singular integrals in Inti.

const ENTITIES

Dictionary mapping EntityKey to GeometricEntity. Contains all entities created in a given session.

SAME_POINTS_TOLERANCE

Two points x and y are considerd the same if norm(x-y) ≤ SAME_POINT_TOLERANCE.

abstract type AbstractDifferentialOperator{N}

A partial differential operator in dimension N.

AbstractDifferentialOperator types are used to define AbstractKernels related to fundamental solutions of differential operators.

abstract type AbstractKernel{T}

A kernel functions K with the signature K(target,source)::T.

See also: SingleLayerKernel, DoubleLayerKernel, AdjointDoubleLayerKernel, HyperSingularKernel

abstract type AbstractMesh{N,T}

An abstract mesh structure in dimension N with primite data of type T (e.g. Float64 for double precision representation).

Concrete subtypes of AbstractMesh should implement ElementIterator for accessing the mesh elements.

See also: Mesh

struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T}

Given an operator Op, construct its free-space adjoint double-layer kernel. This corresponds to the transpose(γ₁,ₓ[G]), where G is the SingleLayerKernel. For operators such as Laplace or Helmholtz, this is simply the normal derivative of the fundamental solution respect to the target variable.

struct BlockArray{T<:StaticArray,N,S} <: AbstractMatrix{T,N}

A struct which behaves like an Array{T,N}, but with the underlying data stored as a Matrix{S}, where S::Number = eltype(T) is the scalar type associated with T. This allows for the use of many blas routines under-the-hood, while providing a convenient interface for handling arrays over StaticArrays.

using StaticArrays
T = SMatrix{2,2,Int,4}
B = Inti.BlockArray{T}([i*j for i in 1:4, j in 1:4])

# output

2×2 Inti.BlockArray{SMatrix{2, 2, Int64, 4}, 2, Int64, 2}:
 [1 2; 2 4]  [3 4; 6 8]
 [3 6; 4 8]  [9 12; 12 16]

struct DimParameters

Parameters associated with the density interpolation method used in bdim_correction.

struct Domain

Representation of a geometrical domain formed by a set of entities with the same geometric dimension. For basic set operations on domains are supported (union, intersection, difference, etc), and they all return a new Domain object.

Calling keys(Ω) returns the set of EntityKeys that make up the domain; given a key, the underlying entities can be accessed with global_get_entity(key).

Domain([f::Function,] keys)

Create a domain from a set of EntityKeys. Optionally, a filter function f can be passed to filter the entities.

Note that all entities in a domain must have the same geometric dimension.

Domain(f::Function, msh::AbstractMesh)

Call Domain(f, ents) on ents = entities(msh).

struct DoubleLayerKernel{T,Op} <: AbstractKernel{T}

Given an operator Op, construct its free-space double-layer kernel. This corresponds to the γ₁ trace of the SingleLayerKernel. For operators such as Laplace or Helmholtz, this is simply the normal derivative of the fundamental solution respect to the source variable.

struct Elastostatic{N,T} <: AbstractDifferentialOperator{N}

Elastostatic operator in N dimensions: -μΔu - (μ+λ)∇(∇⋅u)

Note that the displacement $u$ is a vector of length N since this is a vectorial problem.

struct ElementIterator{E,M} <: AbstractVector{E}

Structure to lazily access elements of type E in a mesh of type M. This is particularly useful for LagrangeElements, where the information to reconstruct the element is stored in the mesh connectivity matrix.

struct EmbeddedQuadrature{L,H,D} <: ReferenceQuadrature{D}

A quadrature rule for the reference shape D based on a high-order quadrature of type H and a low-order quadrature of type L. The low-order quadrature rule is embedded in the sense that its n nodes are exactly the first n nodes of the high-order quadrature rule.

EntityKey

Used to represent the key of a GeometricEntity, comprised of a dim and a tag field, where dim is the geometrical dimension of the entity, and tag is a unique integer identifying the entity.

The sign of the tag field is used to distinguish the orientation of the entity, and is ignored when comparing two EntityKeys for equality.

struct Fejer{N}

N-point Fejer's first quadrature rule for integrating a function over [0,1]. Exactly integrates all polynomials of degree ≤ N-1.

using Inti

q = Inti.Fejer(;order=10)

Inti.integrate(cos,q) ≈ sin(1) - sin(0)

# output

true

struct Gauss{D,N} <: ReferenceQuadrature{D}

Tabulated N-point symmetric Gauss quadrature rule for integration over D.

struct GaussLegendre{N,T}

N-point Gauss-Legendre quadrature rule for integrating a function over [0,1]. Exactly integrates all polynomials of degree ≤ 2N-1.

using Inti

q = Inti.GaussLegendre(;order=10)

Inti.integrate(cos,q) ≈ sin(1) - sin(0)

# output

true

struct GeometricEntity

Geometrical objects such as lines, surfaces, and volumes.

Geometrical entities are stored in a global ENTITIES dictionary mapping EntityKey to the corresponding GeometricEntity, and usually entities are manipulated through their keys.

A GeometricEntity can also contain a pushforward field used to parametrically represent the entry as the image of a reference domain (pushforward.domain) under some function (pushforward.parametrization).

Note that entities are manipulated through their keys, and the GeometricEntity constructor returns the key of the created entity; to retrieve the entity, use the global_get_entity function.

GeometricEntity(shape::String [; translation, rotation, scaling, kwargs...])

Constructs a geometric entity with the specified shape and optional parameters, and returns its key.

Arguments

• shape::String: The shape of the geometric entity.
• translation: The translation vector of the geometric entity. Default is SVector(0, 0, 0).
• rotation: The rotation vector of the geometric entity. Default is SVector(0, 0, 0).
• scaling: The scaling vector of the geometric entity. Default is SVector(1, 1, 1).
• kwargs...: Additional keyword arguments to be passed to the shape constructor.

Supported shapes

• ellipsoid
• torus
• bean
• acorn
• cushion

Helmholtz(; k, dim)

Helmholtz operator in dim dimensions: -Δu - k²u.

The parameter k can be a real or complex number. For purely imaginary wavenumbers, consider using the Yukawa kernel.

struct HyperRectangle{N,T} <: ReferenceInterpolant{ReferenceHyperCube{N},T}

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

struct HyperSingularKernel{T,Op} <: AbstractKernel{T}

Given an operator Op, construct its free-space hypersingular kernel. This corresponds to the transpose(γ₁,ₓγ₁[G]), where G is the SingleLayerKernel. For operators such as Laplace or Helmholtz, this is simply the normal derivative respect to the target variable of the DoubleLayerKernel.

struct IntegralOperator{T} <: AbstractMatrix{T}

A discrete linear integral operator given by

\[I[u](x) = \int_{\Gamma\_s} K(x,y)u(y) ds_y, x \in \Gamma_{t}\]

where $\Gamma_s$ and $\Gamma_t$ are the source and target domains, respectively.

struct IntegralPotential

Represent a potential given by a kernel and a quadrature over which integration is performed.

IntegralPotentials are created using IntegralPotential(kernel, quadrature).

Evaluating an integral potential requires a density σ (defined over the quadrature nodes of the source mesh) and a point x at which to evaluate the integral

\[\int_{\Gamma} K(oldsymbol{x},oldsymbol{y})\sigma(y) ds_y, x \not \in \Gamma\]

Assuming 𝒮 is an integral potential and σ is a vector of values defined on quadrature, calling 𝒮[σ] creates an anonymous function that can be evaluated at any point x.

const LagrangeSquare = LagrangeElement{ReferenceSquare}

struct LagrangeElement{D,Np,T} <: ReferenceInterpolant{D,T}

A polynomial p : D → T uniquely defined by its Np values on the Np reference nodes of D.

The return type T should be a vector space (i.e. support addition and multiplication by scalars). For istance, T could be a number or a vector, but not a Tuple.

const LagrangeLine = LagrangeElement{ReferenceLine}

const LagrangeSquare = LagrangeElement{ReferenceSquare}

const LagrangeTetrahedron = LagrangeElement{ReferenceTetrahedron}

const LagrangeTriangle = LagrangeElement{ReferenceTriangle}

Laplace(; dim)

Laplace's differential operator in dim dimension: $-Δu$. ```

Note the negative sign in the definition.

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

Unstructured mesh defined by a set of nodes(of typeSVector{N,T}`), and a dictionary mapping element types to connectivity matrices. Each columns of a given connectivity matrix stores the integer tags of the nodes in the mesh comprising the element.

Additionally, the mesh contains a mapping from EntityKeys to the tags of the elements composing the entity. This can be used to extract submeshes from a given mesh using e.g. view(msh,Γ) or msh[Γ], where Γ is a Domain.

See elements for a way to iterate over the elements of a mesh.

const ModifiedHelmholtz

Type alias for the Yukawa operator.

MultiIndex{N}

Wrapper around NTuple{N,Int} mimicking a multi-index in ℤ₀ᴺ.

ParametricElement{D,T,F} <: ReferenceInterpolant{D,T}

An element represented through a explicit function f mapping D into the element. For performance reasons, f should take as input a StaticVector and return a StaticVector or StaticArray.

See also: ReferenceInterpolant, LagrangeElement

ParametricElement(f, d::HyperRectangle)

Construct the element defined as the image of f over d.

struct PolynomialSpace{D,K}

The space of all polynomials of degree ≤K, commonly referred to as ℙₖ.

The type parameter D, of singleton type, is used to determine the reference domain of the polynomial basis. In particular, when D is a hypercube in d dimensions, the precise definition is ℙₖ = span{𝐱ᶿ : 0≤max(θ)≤ K}; when D is a d-dimensional simplex, the space is ℙₖ = span{𝐱ᶿ : 0≤sum(θ)≤ K}, where θ ∈ 𝐍ᵈ is a multi-index.

See also: monomial_basis, lagrange_basis

struct Quadrature{N,T} <: AbstractVector{QuadratureNode{N,T}}

A collection of QuadratureNodes used to integrate over an AbstractMesh.

Quadrature(Ω::Domain; meshsize, qorder)

Construct a Quadrature over the domain Ω with a mesh of size meshsize and quadrature order qorder.

Quadrature(msh::AbstractMesh, etype2qrule::Dict)
Quadrature(msh::AbstractMesh, qrule::ReferenceQuadrature)
Quadrature(msh::AbstractMesh; qorder)

Construct a Quadrature for msh, where for each element type E in msh the reference quadrature q = etype2qrule[E] is used. When a single qrule is passed, it is used for all element types in msh.

If an order keyword is passed, a default quadrature of the desired order is used for each element type usig _qrule_for_reference_shape.

For co-dimension one elements, the normal vector is also computed and stored in the QuadratureNodes.

QuadratureNode{N,T<:Real}

A point in ℝᴺ with a weight for performing numerical integration. A QuadratureNode can optionally store a normal vector.

const ReferenceCube = ReferenceHyperCube{3}

Singleton type representing the unit cube [0,1]³.

struct ReferenceHyperCube{N} <: ReferenceShape{N}

Singleton type representing the axis-aligned hypercube in N dimensions with the lower corner at the origin and the upper corner at (1,1,…,1).

abstract type ReferenceInterpolant{D,T}

Interpolanting function mapping points on the domain D<:ReferenceShape (of singleton type) to a value of type T.

Instances el of ReferenceInterpolant are expected to implement:

• el(x̂): evaluate the interpolation scheme at the (reference) coordinate x̂ ∈ D.
• jacobian(el,x̂) : evaluate the jacobian matrix of the interpolation at the (reference) coordinate x ∈ D.

const ReferenceLine = ReferenceHyperCube{1}

Singleton type representing the [0,1] segment.

abstract type ReferenceQuadrature{D}

A quadrature rule for integrating a function over the domain D <: ReferenceShape.

Calling x,w = q() returns the nodes x, given as SVectors, and weights w, for performing integration over domain(q).

abstract type ReferenceShape

A fixed reference domain/shape. Used mostly for defining more complex shapes as transformations mapping an ReferenceShape to some region of ℜᴹ.

See e.g. ReferenceLine or ReferenceTriangle for some examples of concrete subtypes.

struct ReferenceSimplex{N}

Singleton type representing the N-simplex with N+1 vertices (0,...,0),(0,...,0,1),(0,...,0,1,0),(1,0,...,0)

const ReferenceSquare = ReferenceHyperCube{2}

Singleton type representing the unit square [0,1]².

struct ReferenceTetrahedron

Singleton type representing the tetrahedron with vertices (0,0,0),(0,0,1),(0,1,0),(1,0,0)

struct ReferenceTriangle

Singleton type representing the triangle with vertices (0,0),(1,0),(0,1)

struct SingleLayerKernel{T,Op} <: AbstractKernel{T}

The free-space single-layer kernel (i.e. the fundamental solution) of an Op <: AbstractDifferentialOperator.

Stokes(; μ, dim)

Stokes operator in dim dimensions: $[-μΔu + ∇p, ∇⋅u]$.

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

View into a parent mesh over a given domain.

A submesh implements the interface for AbstractMesh; therefore you can iterate over elements of the submesh just like you would with a mesh.

Construct SubMeshs using view(parent,Ω::Domain).

TensorProductQuadrature{N,Q}

A tensor-product of one-dimension quadrature rules. Integrates over [0,1]^N.

Examples

qx = Inti.Fejer(10)
qy = Inti.Fejer(15)
q  = Inti.TensorProductQuadrature(qx,qy)

struct VioreanuRokhlin{D,N} <: ReferenceQuadrature{D}

Tabulated N-point Vioreanu-Rokhlin quadrature rule for integration over D.

Yukawa(; λ, dim)

Yukawa operator, also known as modified Helmholtz, in dim dimensions: $-Δu + λ²u$.

The parameter λ is a positive number. Note the negative sign in front of the Laplacian.

iterate(Ω::Domain)

Iterating over a domain means iterating over its entities.

DoubleLayerPotential(op::AbstractDifferentialOperator, source::Quadrature)

An IntegralPotential over source with kernel given by DoubleLayerKernel(op).

SingleLayerPotential(op::AbstractDifferentialOperator, source::Quadrature)

An IntegralPotential over source with kernel given by SingleLayerKernel(op).

_blocksize_normalized(A::BlockArray)

Like blocksize, but appends 1s if A is a higher-dimensional.

For example, a BlockArray{SVector{3,Float64}, 2} has a blocksize of (3,), but a normalized_blocksize of (3, 1).

_get_gauss_and_qweights(R::Type{<:ReferenceShape{D}}, N) where D

Returns the N-point symmetric gaussian qnodes and qweights (x, w) for integration over R.

_get_vioreanurokhlin_qcoords_and_qweights(R::Type{<:ReferenceShape{D}}, N) where D

Returns the N-point Vioreanu-Rokhlin qnodes and qweights (x, w) for integration over R.

_green_multiplier(s::Symbol)

Return -1.0 if s == :inside, 0.0 if s == :outside, and -0.5 if s == :on; otherwise, throw an error. The orientation is relative to the normal of the bounding curve/surface.

_green_multiplier(x, quad)

Helper function to help determine the constant σ in the Green identity S[γ₁u](x)

• D[γ₀u](x) + σ*u(x) = 0. This can be used as a predicate to determine whether a

point is inside a domain or not.

_meshgen(f,d::HyperRectangle,sz)

Create prod(sz) elements of ParametricElement type representing the push forward of f on each of the subdomains defined by a uniform cartesian mesh of d of size sz.

_normal(jac::SMatrix{M,N}, s = 1)

Given a an M by N matrix representing the jacobian of a codimension one object, compute the normal vector. If s=-1, the normal vector is flipped.

_qrule_for_reference_shape(ref,order)

Given a reference shape and a desired quadrature order, return an appropiate quadrature rule.

acorn(; translation, rotation, scaling, labels)

Create an acorn entity in 3D, and apply optional transformations. Returns the key.

adaptive_correction(iop::IntegralOperator; [maxdist, rtol, threads = true, kwargs...])
adaptive_correction(iop::IntegralOperator, maxdist, quads_dict::Dict, threads = true)

This function computes a sparse correction for the integral operator iop, addressing its singular or nearly singular entries.

The parameter maxdist specifies the maximum distance between target points and source elements to be considered for correction (only interactions within this distance are corrected).

The parameter rtol defines the tolerance for the adaptive quadrature used to compute the corrections for singular or nearly singular entries.

Additional kwargs arguments are passed to adaptive_quadrature; see its documentation for more information.

Selecting maxdist and rtol involves balancing accuracy and computational cost. A smaller maxdist reduces the number of corrections but may impact accuracy. Conversely, a smaller rtol improves correction accuracy but increases computational expense. The ideal values for maxdist and rtol depend on the kernel and the mesh/quadrature rule applied.

By default, maxdist and rtol are estimated using the local_correction_dist_and_tol, but it is often possible to improve performance by manually tunning these parameters.

Advanced usage

For finer control, you can provide a dictionary quads_dict that contains quadrature rules for each reference element type present in the mesh of source(iop). This allows you to fine-tune the quadrature rules for specific element types (e.g. use a fixed quadrature rule instead of an adaptive one).

The dictionary quads_dict must adhere to the following structure:

• quads_dict[E].nearfield_quad: A function that integrates over the nearfield of the reference element type E. Used in the nearly-singular correction.
• quads_dict[E].radial_quad: A function that integrates over the radial direction of the reference element type E. Used in the singular correction.
• quads_dict[E].angular_quad: A function that integrates over the angular direction of the reference element type E. Used in the singular correction.

Here is an example of how to implement a custom quads_dict given an iop:

quads_dict = Dict()
msh = Inti.mesh(source(iop))
for E in Inti.element_types(msh)
    ref_domain = Inti.reference_domain(E)
    quads = (
        nearfield_quad = Inti.adaptive_quadrature(ref_domain; atol),
        radial_quad    = Inti.GaussLegendre(;order=5),
        angular_quad   = Inti.GuassLegendre(;order=20),
    )
    quads_dict[E] = quads
end

This will use an adaptive quadrature rule for the nearfield and fixed Gauss-Legendre quadrature rules for the radial and angular directions when computing the singular correction in polar coordinates on the reference domain. You can then call adaptive_correction(iop, maxdist, quads_dict) to use the custom quadrature.

adaptive_quadrature(ref_domain::ReferenceShape; kwargs...)

Return a function quad callable as quad(f) that integrates the function f over the reference shape ref_domain. The keyword arguments are passed to HAdaptiveIntegration.integrate.

adj_double_layer_hypersingular(; op, target, source, compression,
correction)

Similar to single_double_layer, but for the adjoint double-layer and hypersingular operators. See the documentation of [single_double_layer] for a description of the arguments.

ambient_dimension(x)

Dimension of the ambient space where x lives. For geometrical objects this can differ from its geometric_dimension; for example a triangle in ℝ³ has ambient dimension 3 but geometric dimension 2, while a curve in ℝ³ has ambient dimension 3 but geometric dimension 1.

assemble_fmm(iop; rtol)

Set up a 2D or 3D FMM for evaluating the discretized integral operator iop associated with the op. In 2D the FMM2D or FMMLIB2D library is used (whichever was most recently loaded) while in 3D FMM3D is used.

assemble_hmatrix(iop[; atol, rank, rtol, eta])

Assemble an H-matrix representation of the discretized integral operator iop using the HMatrices.jl library.

See the documentation of HMatrices for more details on usage and other keyword arguments.

assemble_matrix(iop::IntegralOperator; threads = true)

Assemble a dense matrix representation of an IntegralOperator.

bdim_correction(op,X,Y,S,D; green_multiplier, kwargs...)

Given a op and a (possibly inaccurate) discretizations of its single and double-layer operators S and D (taking a vector of values on Y and returning a vector on of values on X), compute corrections δS and δD such that S + δS and D + δD are more accurate approximations of the underlying single- and double-layer integral operators.

See [8] for more details on the method.

Arguments

Required:

• op must be an AbstractDifferentialOperator
• Y must be a Quadrature object of a closed surface
• X is either inside, outside, or on Y
• S and D are approximations to the single- and double-layer operators for op taking densities in Y and returning densities in X.
• green_multiplier (keyword argument) is a vector with the same length as X storing the value of μ(x) for x ∈ X in the Green identity S\[γ₁u\](x) - D\[γ₀u\](x) + μ*u(x) = 0. See _green_multiplier.

Optional kwargs:

• parameters::DimParameters: parameters associated with the density interpolation method
• derivative: if true, compute the correction to the adjoint double-layer and hypersingular operators instead. In this case, S and D should be replaced by a (possibly innacurate) discretization of adjoint double-layer and hypersingular operators, respectively.
• maxdist: distance beyond which interactions are considered sufficiently far so that no correction is needed. This is used to determine a threshold for nearly-singular corrections when X and Y are different surfaces. When X === Y, this is not needed.

bean(; translation, rotation, scaling, labels)

Create a bean entity in 3D, and apply optional transformations. Returns the key.

blocksize(A::BlockArray)

The size of an individual entry of A.

boundary(Ω::Domain)

Return the external boundaries of a domain.

See also: external_boundary, internal_boundary, skeleton.

boundary_idxs(el::LagrangeElement)

The indices of the nodes in el that define the boundary of the element.

cart2sph(x,y,z)

Map cartesian coordinates x,y,z to spherical ones r, θ, φ representing the radius, elevation, and azimuthal angle respectively. The convention followed is that 0 ≤ θ ≤ π and -π < φ ≤ π. Same as the cart2sph function in MATLAB.

connectivity(msh::AbstractMesh,E::DataType)

Return the connectivity matrix for elements of type E in msh. The integer tags in the matrix refer to the points in nodes(msh)

coords(q)

Return the spatial coordinates of q.

cushion(; translation, rotation, scaling, labels)

Create a cushion entity in 3D, and apply optional transformations. Returns the key.

degree(el::LagrangeElement)
degree(el::Type{<:LagrangeElement})

The polynomial degree el.

dimension(space)

The length of a basis for space; i.e. the number of linearly independent elements required to span space.

dom2elt(m::Mesh,Ω,E)::Vector{Int}

Compute the element indices idxs of the elements of type E composing Ω.

dom2qtags(Q::Quadrature, dom::Domain)

Given a domain, return the indices of the quadratures nodes in Q associated to its quadrature.

domain(f)

Given a function-like object f: Ω → R, return Ω.

domain(msh::AbstractMesh)

Return a [Domain] containing of all entities covered by the mesh.

domain(Q::Quadrature)

The Domain over which Q performs integration.

domain(q::ReferenceQuadrature)

The domain of integratino for quadrature rule q.

element_types(msh::AbstractMesh)

Return the element types present in the msh.

elements(msh::AbstractMesh [, E::DataType])

Return the elements of a msh. Passing and element type E will restricts to elements of that type.

A common pattern to avoid type-instabilies in performance critical parts of the code is to use a function barrier, as illustrated below:

for E in element_types(msh)
    _long_computation(elements(msh, E), args...)
end

@noinline function _long_computation(iter, args...)
    for el in iter # the type of el is known at compile time
        # do something with el
    end
end

where a dynamic dispatch is performed only on the element types (typically small for a given mesh).

ellipsoid(; translation, rotation, scaling, labels)

Create an ellipsoid entity in 3D, and apply optional transformations. Returns the key of the created entity.

ent2etags(msh::AbstractMesh)

Return a dictionary mapping entities to a dictionary of element types to element tags.

entities(Ω::Domain)

Return all entities making up a domain (as a set of EntityKeys).

etype_to_nearest_points(X,Y::Quadrature; maxdist)

For each element el in Y.mesh, return a list with the indices of all points in X for which el is the nearest element. Ignore indices for which the distance exceeds maxdist.

external_boundary(Ω::Domain)

Return the external boundaries inside a domain. These are entities in the skeleton of Ω which are not in the internal boundaries of Ω.

See also: internal_boundary, skeleton.

fibonnaci_points_sphere(N,r,c)

Return N points distributed (roughly) in a uniform manner on the sphere of radius r centered at c.

flip_normal(q::QuadratureNode)

Return a new QuadratureNode with the normal vector flipped.

gauss_curvature(Q::Quadrature)

Compute the gauss_curvature at each quadrature node in Q.

gauss_curvature(τ, x̂)

Calculate the Gaussian curvature of the element τ at the parametric coordinate x̂.

geometric_dimension(x)

NNumber of degrees of freedom necessary to locally represent the geometrical object. For example, lines have geometric dimension of 1 (whether in ℝ² or in ℝ³), while surfaces have geometric dimension of 2.

global_get_entity(k::EntityKey)

Retrieve the GeometricEntity corresponding to the EntityKey k from the global ENTITIES dictionary.

gmsh_curve(f::Function, a, b; npts=100, meshsize = 0, tag=-1)

Create a curve in the current gmsh model given by {f(t) : t ∈ (a,b) } where f is a function from ℝ to ℝ^3. The curve is approximated by C² b-splines passing through npts equispaced in parameter space. If a meshsize is given, gmsh will use it when meshing the curve.

hesssian(el,x)

Given a (possibly vector-valued) functor f : 𝐑ᵐ → 𝐅ⁿ, return the n × m × m matrix Aᵢⱼⱼ = ∂²fᵢ/∂xⱼ∂xⱼ. By default ForwardDiff is used to compute the hessian, but you should overload this method for specific f if better performance and/or precision is required.

Note: both x and f(x) are expected to be of SVector type.

image(f)

Given a function-like object f: Ω → R, return f(Ω).

import_mesh(filename = nothing; dim=3)

Open filename and create a Mesh from the gmsh model in it.

If filename is nothing, the current gmsh model is used. Note that this assumes that the Gmsh API has been initialized through gmsh.initialize.

Passing dim=2 will create a two-dimensional mesh by projecting the original mesh onto the x,y plane.

integrate(f,quad::Quadrature)

Compute ∑ᵢ f(qᵢ)wᵢ, where the qᵢ are the quadrature nodes of quad, and wᵢ are the quadrature weights.

Note that you must define f(::QuadratureNode): use q.coords and q.normal if you need to access the coordinate or normal vector at que quadrature node.

integrate(f,q::ReferenceQuadrature)
integrate(f,x,w)

Integrate the function f using the quadrature rule q. This is simply sum(f.(x) .* w), where x and w are the quadrature nodes and weights, respectively.

The function f should take an SVector as input.

integrate_with_error_estimate(f, quad::EmbeddedQuadrature, norm = LinearAlgebra.norm)

Return I, E where I is the estimated integral of f over domain(quad) using the high-order quadrature and E is the error estimate obtained by taking the norm of the difference between the high and low-order quadratures in quad.

integration_measure(f, x̂)

Given the Jacobian matrix J of a transformation f : ℝᴹ → ℝᴺ compute the integration measure √det(JᵀJ) at the parametric coordinate x̂

interface_method(x)

A method of an abstract type for which concrete subtypes are expected to provide an implementation.

internal_boundary(Ω::Domain)

Return the internal boundaries of a Domain. These are entities in skeleton(Ω) which appear at least twice as a boundary of entities in Ω.

interpolation_order(qrule::ReferenceQuadrature)

The interpolation order of a quadrature rule is defined as the the smallest k such that there exists a unique polynomial in PolynomialSpace{D,k} that minimizes the error in approximating the function f at the quadrature nodes.

For an N-point Gauss quadrature rule on the segment, the interpolation order is N-1 since N points uniquely determine a polynomial of degree N-1.

For a triangular reference domain, the interpolation order is more difficult to define. An unisolvent three-node quadrature on the triangular, for example, has an interpolation order k=1 since the three nodes uniquely determine a linear polynomial, but a four-node quadrature may also have an interpolation order k=1 since for k=2 there are multiple polynomials that pass through the four nodes.

jacobian(f,x)

Given a (possibly vector-valued) functor f : 𝐑ᵐ → 𝐅ⁿ, return the n × m matrix Aᵢⱼ = ∂fᵢ/∂xⱼ. By default ForwardDiff is used to compute the jacobian, but you should overload this method for specific f if better performance and/or precision is required.

Note: both x and f(x) are expected to be of SVector type.

kress_change_of_variables(P)

Return a change of variables mapping [0,1] to [0,1] with the property that the first P-1 derivatives of the transformation vanish at x=0.

kress_change_of_variables_periodic(P)

Like kress_change_of_variables, this change of variables maps the interval [0,1] onto itself, but the first P derivatives of the transformation vanish at both endpoints (thus making it a periodic function).

This change of variables can be used to periodize integrals over the interval [0,1] by mapping the integrand into a new integrand that vanishes (to order P) at both endpoints.

lagrange_basis(nodes,[sp::AbstractPolynomialSpace])

Return the set of n polynomials in sp taking the value of 1 on node i and 0 on nodes j ≂̸ i for 1 ≤ i ≤ n.

lagrange_basis(qrule::ReferenceQuadrature)

Return a function L : ℝᴺ → ℝᵖ where N is the dimension of the domain of qrule, and p is the number of nodes in qrule. The function L is a polynomial in polynomial_space(qrule), and L(xⱼ)[i] = δᵢⱼ (i.e. the ith component of L is the ith Lagrange basis).

lagrange_basis(E::Type{<:LagrangeElement})

Return the Lagrange basis B for the element E. Evaluating B(x) yields the value of each basis function at x.

laurent_coefficients(f, h, order) --> f₋₂, f₋₁, f₀

Given a one-dimensional function f, return f₋₂, f₋₁, f₀ such that f(x) = f₋₂ / x^2 + f₋₁ / x + f₀ + 𝒪(x) as x -> 0, where we assume that fₙ = 0 for n < N.

The order argument is an integer that indicates the order of the singularity at the origin:

• Val{-2}: The function has a singularity of order -2 at the origin.
• Val{-1}: The function has a singularity of order -1 at the origin, so f₋₂ = 0.
• Val{0}: The function has a finite part at the origin, so f₋₂ = f₋₁ = 0.

line(a,b)

Create a [GeometricEntity] representing a straight line connecting points a and b. The points a and b can be either SVectors or a Tuple.

The parametrization of the line is given by f(u) = a + u(b - a), where 0 ≤ u ≤ 1.

local_correction_dist_and_tol(iop::IntegralOperator, kmax = 10, ratio = 8)

Try to estimate resonable maxdist and rtol parameters for the adaptive_correction function, where maxdist is at most kmax times the radius of the largest element in the source mesh of iop. See the Extended help for more details.

Extended help

The heuristic works as follows, where we let K = kernel(iop) and msh = mesh(source(iop)):

1. Pick the largest element in msh
2. Let h be the radius of el
3. For k between 1 and kmax, estimate the (relative) quadrature error when integrating y -> K(x,y) for x at a distance k * h from the center of the element using a regular quadrature rule
4. Find a k such that ratio between errors at distances k * h and (k + 1) * h is below ratio. This indicates stagnation in the error, and suggests that little is gained by increasing the distance.
5. Return maxdist = k * h and rtol as the error at distance k * h.

mean_curvature(Q::Quadrature)

Compute the mean_curvature at each quadrature node in Q.

mean_curvature(τ, x̂)

Calculate the mean curvature of the element τ at the parametric coordinate x̂.

measure(k::EntityKey, rtol)

Compute the length/area/volume of the entity k using an adaptive quadrature with a relative tolerance rtol. Assumes that the entity has an explicit parametrization.

meshgen!(mesh,Ω,sz)

Similar to meshgen, but append entries to mesh.

meshgen(Ω, n; T = Float64)
meshgen(Ω, n_dict; T = Float64)
meshgen(Ω; meshsize, T = Float64)

Generate a Mesh for the domain Ω where each curve is meshed using n elements. Passing a dictionary allows for a finer control; in such cases, n_dict[ent] should return an integer for each entity ent in Ω of geometric_dimension one.

Alternatively, a meshsize can be passed, in which case, the number of elements is computed as so as to obtain an average mesh size of meshsize. Note that the actual mesh size may vary significantly for each element if the parametrization is far from uniform.

The mesh is created with primitive data of type T.

This function requires the entities forming Ω to have an explicit parametrization.

monomial_basis(sp::PolynomialSpace)

Return a function f : ℝᴺ → ℝᵈ, where N is the dimension of the domain of sp containing a basis of monomials 𝐱ᶿ spanning the polynomial space PolynomialSpace.

near_interaction_list(X,Y::AbstractMesh; tol)

For each element el of type E in Y, return the indices of the points in X which are closer than tol to the center of el.

This function returns a dictionary where e.g. dict[E][5] --> Vector{Int} gives the indices of points in X which are closer than tol to the center of the fifth element of type E.

If tol is a Dict, then tol[E] is the tolerance for elements of type E.

new_tag(dim)

Return a new tag for an entity of dimension dim so that EntityKey(dim, tag) is not already in ENTITIES.

nodes(msh::SubMesh)

A view of the nodes of the parent mesh belonging to the submesh. The ordering is given by the nodetags function.

nodetags(msh::SubMesh)

Return the tags of the nodes in the parent mesh belonging to the submesh.

normal(el, x̂)

Return the normal vector of el at the parametric coordinate x̂.

normal(q)

Return the normal vector of q, if it exists.

notimplemented()

Things which should probably be implemented at some point.

order(q::ReferenceQuadrature)

A quadrature of order p (sometimes called degree of precision) integrates all polynomials of degree ≤ p but not ≤ p + 1.

order(el::LagrangeElement)

The order of the element's interpolating polynomial (e.g. a LagrangeLine with 2 nodes defines a linear polynomial, and thus has order 1).

orientation(msh::AbstractMesh,E::DataType)

Return the orientation of the elements of type E in msh (1 if normal and -1 if inverted).

parametric_curve(f, a::Real, b::Real)

Create a [GeometricEntity] representing a parametric curve defined by the {f(t) | a ≤ t ≤ b}. The function f should map a scalar to an SVector.

Flipping the orientation is supported by passing a > b.

    parametric_surface(f, lc, hc, boundary = nothing; kwargs...)

Create a parametric surface defined by the function f over the rectangular domain defined by the lower corner lc and the upper corner hc. The optional boundary argument can be used to specify the boundary curves of the surface.

Arguments

• f: A function that takes two arguments x and y and returns a tuple (u, v) representing the parametric coordinates of the surface at (x, y).
• lc: A 2-element array representing the lower corner of the rectangular domain.
• hc: A 2-element array representing the upper corner of the rectangular domain.
• boundary: An optional array of boundary curves that define the surface.

Keyword Arguments

• kwargs: Additional keyword arguments that can be passed to the GeometricEntity constructor.

Returns

• The key of the created GeometricEntity.

polar_decomposition(shape::ReferenceSquare, x̂::SVector{2,Float64})

Decompose the square [0,1] × [0,1] into four triangles, and return four tuples of the form θₛ, θₑ, ρ where θₛ and θₑ are the initial and final angles of the triangle, and ρ is the function that gives the distance from x̂ to the border of the square in the direction θ.

polar_decomposition(shape::ReferenceTriangle, x̂::SVector{2,Float64})

Decompose the triangle {x,y ≥ 0, x + y ≤ 1} into three triangles, and return three tuples of the form θₛ, θₑ, ρ where θₛ and θₑ are the initial and final angles of the triangle, and ρ is the function that gives the distance from x̂ to the border of the triangle in the direction θ.

polynomial_solutions_vdim(op, order[, center])

For every monomial term pₙ of degree order, compute a polynomial Pₙ such that ℒ[Pₙ] = pₙ, where ℒ is the differential operator associated with op. This function returns {pₙ,Pₙ,γ₁Pₙ}, where γ₁Pₙ is the generalized Neumann trace of Pₙ.

Passing a point center will shift the monomials and solutions accordingly.

polynomial_space(qrule::ReferenceQuadrature)

Return a PolynomialSpace associated with the interpolation_order of the quadrature nodes of qrule.

qcoords(q)

Return the coordinate of the quadrature nodes associated with q.

quadrature_to_node_vals(Q::Quadrature, qvals::AbstractVector)

Given a vector qvals of scalar values at the quadrature nodes of Q, return a vector ivals of scalar values at the interpolation nodes of Q.mesh.

qweights(q)

Return the quadrature weights associated with q.

reference_nodes(el::LagrangeElement)
reference_nodes(::Type{<:LagrangeElement})

Return the reference nodes on domain(el) used for the polynomial interpolation. The function values on these nodes completely determines the interpolating polynomial.

We use the same convention as gmsh for defining the reference nodes and their order (see node ordering on gmsh documentation).

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 knonw return type should extend return_type(::T,args...) to avoid relying on promote_op.

rotation_matrix(rot)

Constructs a rotation matrix given the rotation angles around the x, y, and z axes.

Arguments

• rot: A tuple or vector containing the rotation angles in radians for each axis.

Returns

• R::SMatrix: The resulting rotation matrix.

single_double_layer(; op, target, source::Quadrature, compression,
correction, derivative = false)

Construct a discrete approximation to the single- and double-layer integral operators for op, mapping values defined on the quadrature nodes of source to values defined on the nodes of target. If derivative = true, return instead the adjoint double-layer and hypersingular operators (which are the generalized Neumann trace of the single- and double-layer, respectively).

For finer control, you must choose a compression method and a correction method, as described below.

Compression

The compression argument is a named tuple with a method field followed by method-specific fields. It specifies how the dense linear operators should be compressed. The available options are:

• (method = :none, ): no compression is performed, the resulting matrices are dense. This is the default, but not recommended for large problems.
• (method =:hmatrix, tol): the resulting operators are compressed using hierarchical matrices with an absolute tolerance tol (defaults to 1e-8).
• (method = :fmm, tol): the resulting operators are compressed using the fast multipole method with an absolute tolerance tol (defaults to 1e-8).

Correction

The correction argument is a named tuple with a method field followed by method-specific fields. It specifies how the singular and nearly-singular integrals should be computed. The available options are:

• (method = :none, ): no correction is performed. This is not recommended, as the resulting approximation will be inaccurate if the kernel is singular and source and target are not sufficiently far from each other.
• (method = :adaptive, maxdist, tol): correct interactions corresponding to entries of target and elements of source that are within maxdist of each other. The singular (including finite part) interactions are computed in polar coordinates, while the near-singular interactions are computing using an adaptive quadrature rule. The tol argument specifies the tolerance of the adaptive integration. See adaptive_correction for more details.
• (method = :dim, maxdist, target_location): use the density interpolation method to compute the correction. maxdist specifies the distance between source and target points above which no correction is performed (defaults to Inf). target_location should be either :inside, :outside, or :on, and specifies where the targetpoints lie relative to the to thesourcecurve/surface (which is assumed to be closed). Whentarget === source,targetlocationis not needed. See [bdimcorrection](@ref) and [vdim_correction`] for more details.

single_double_layer_potential(; op, source)

Return the single- and double-layer potentials for op as IntegralPotentials.

singularity_order(K)

Given a kernel K with signature K(target,source)::T, return the order of the singularity of K at target = source. Order n means that K(x,y) ∼ (x - y)^n as x -> y.

skeleton(Ω::Domain)

Return all the boundaries of the domain, i.e. the domain's skeleton.

stack_weakdeps_env!(; verbose = false, update = false)

Push to the load stack an environment providing the weak dependencies of Inti.jl. This allows benefiting from additional functionalities of Inti.jl which are powered by weak dependencies without having to manually install them in your environment.

Set update=true if you want to update the weakdeps environment.

Examples:

Inti.stack_weakdeps_env!()
using HMatrices

standard_basis_vector(k, ::Val{N})

Create an SVector of length N with a 1 in the kth position and zeros elsewhere.

svector(f,n)

Create an SVector of length n, computing each element as f(i), where i is the index of the element.

torus(; r, R, translation, rotation, scaling, labels)

Create a torus entity in 3D, and apply optional transformations. Returns the key. The parameters r and R are the inner and outer radii of the torus.

uniform_points_circle(N,r,c)

Return N points uniformly distributed on a circle of radius r centered at c.

vdim_correction(op,X,Y,Y_boundary,S,D,V; green_multiplier, kwargs...)

Compute a correction to the volume potential V : Y → X such that V + δV is a more accurate approximation of the underlying volume potential operator. The correction is computed using the (volume) density interpolation method.

This function requires a op::AbstractDifferentialOperator, a target set X, a source quadrature Y, a boundary quadrature Y_boundary, approximations S : Y_boundary -> X and D : Y_boundary -> X to the single- and double-layer potentials (correctly handling nearly-singular integrals), and a naive approximation of the volume potential V. The green_multiplier is a vector of the same length as X storing the value of μ(x) for x ∈ X in the Green identity (see _green_multiplier).

See [11] for more details on the method.

Optional kwargs:

• interpolation_order: the order of the polynomial interpolation. By default, the maximum order of the quadrature rules is used.
• maxdist: distance beyond which interactions are considered sufficiently far so that no correction is needed. This is used to determine a threshold for nearly-singular corrections.
• center: the center of the basis functions. By default, the basis functions are centered at the origin.
• shift: a boolean indicating whether the basis functions should be shifted and rescaled to each element.

vdim_mesh_center(msh)

Point x which minimizes ∑ (x-xⱼ)²/r²ⱼ, where xⱼ and rⱼ are the circumcenter and circumradius of the elements of msh, respectively.

vertices(el::LagrangeElement)

Coordinates of the vertices of el.

vertices_idxs(el::LagrangeElement)

The indices of the nodes in el that define the vertices of the element.

volume_potential(; op, target, source::Quadrature, compression, correction)

Compute the volume potential operator for a given PDE.

Arguments

• op: The PDE (Partial Differential Equation) to solve.
• target: The target domain where the potential is computed.
• source: The source domain where the potential is generated.
• compression: The compression method to use for the potential operator.
• correction: The correction method to use for the potential operator.

Returns

The volume potential operator V that represents the interaction between the target and source domains.

Compression

The compression argument is a named tuple with a method field followed by method-specific fields. It specifies how the dense linear operators should be compressed. The available options are:

• (method = :none, ): no compression is performed, the resulting matrices are dense.
• (method =:hmatrix, tol): the resulting operators are compressed using hierarchical matrices with an absolute tolerance tol (defaults to 1e-8).
• (method = :fmm, tol): the resulting operators are compressed using the fast multipole method with an absolute tolerance tol (defaults to 1e-8).

Correction

The correction argument is a named tuple with a method field followed by method-specific fields. It specifies how the singular and nearly-singular integrals should be computed. The available options are:

• (method = :none, ): no correction is performed. This is not recommented, as the resulting approximation will be inaccurate if the source and target are not sufficiently far apart.
• (method = :dim, maxdist, target_location): use the density interpolation method to compute the correction. maxdist specifies the distance between source and target points above which no correction is performed (defaults to Inf). target_location should be either :inside, :outside, or :on, and specifies where the targetpoints lie relative to the to thesource's boundary. Whentarget === source,target_location` is not needed.

Details

The volume potential operator is computed by assembling the integral operator V using the single-layer kernel G. The operator V is then compressed using the specified compression method. If no compression is specified, the operator is returned as is. If a correction method is specified, the correction is computed and added to the compressed operator.