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, :hcubature]

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 blas routines under-the-hood, while providing a convenient interface for handling matrices 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}:
 [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.

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

N-point Kronrod rule obtained by adding n+1 points to a Gauss quadrature containing n points. The order is either 3n + 1 for n even or 3n + 2 for n odd.

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.

_copyto!(target,source)

Defaults to Base.copyto!, but includes some specialized methods to copy from a Matrix of SMatrix to a Matrix of Numbers and viceversa.

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

Given a an M by N matrix representing the jacobian of a codimension one object, compute the normal vector.

_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; tol, maxdist = farfield_distance(iop; atol), maxsplit = 1000])

Given an integral operator iop, this function provides a sparse correction to iop for the entries i,j such that the distance between the i-th target and the j-th source is less than maxdist.

Choosing maxdist is a trade-off between accuracy and efficiency. The smaller the value, the fewer corrections are needed, but this may compromise the accuracy. For a fixed quadrature, the size of maxdist has to grow as the tolerance tol decreases. The default [farfield_distance(iop; tol)](@ref) provides a heuristic to determine a suitablemaxdist`.

The correction is computed by using the adaptive_integration routine, with a tolerance atol and a maximum number of subdivisions maxsplit; see adaptive_integration for more details.

adaptive_integration(f, τ̂::RefernceShape; kwargs...)
adaptive_integration(f, qrule::EmbeddedQuadrature; kwargs...)

Use an adaptive procedure to estimate the integral of f over τ̂ = domain(qrule). The following optional keyword arguments are available:

• atol::Real=0.0: absolute tolerance for the integral estimate
• rtol::Real=0.0: relative tolerance for the integral estimate
• maxsplit::Int=1000: maximum number of times to split the domain
• norm::Function=LinearAlgebra.norm: norm to use for error estimates
• buffer::BinaryHeap: a pre-allocated buffer to use for the adaptive procedure (see allocate_buffer)

adaptive_integration_singular(f, τ̂, x̂ₛ; kwargs...)

Similar to adaptive_integration, but indicates that f has an isolated (integrable) singularity at x̂ₛ ∈ x̂ₛ.

The integration is performed by splitting τ̂ so that x̂ₛ is a fixed vertex, guaranteeing that f is never evaluated at x̂ₛ. Aditionally, a suitable change of variables may be applied to alleviate the singularity and improve the rate of convergence.

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.

allocate_buffer(f, quad::EmbeddedQuadrature)

Create the buffer needed for the call adaptive_integration(f, τ̂; buffer, ...).

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; atol)

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 innacurate) 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 [7] 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.

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.

decompose(s::ReferenceShape,x)

Decompose an ReferenceShape into LagrangeElements so that x is a fixed vertex of the children elements.

The decomposed elements may be oriented differently than the parent, and thus care has to be taken regarding e.g. normal vectors.

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.

farfield_distance(iop::IntegralOperator; tol, maxiter = 10)
farfield_distance(K, Q::Quadrature; tol, maxiter = 10)

Return an estimate of the distance d such that the (absolute) quadrature error of the integrand y -> K(x,y) is below tol for x at a distance d from the center of the largest element in Q; when an integral operator is passed, we have Q::Quadrature = source(iop) and K = kernel(iop).

The estimate is computed by finding the first integer n such that the quadrature error on the largest element τ lies below tol for points x satisfying dist(x,center(τ)) = n*radius(τ).

Note that the desired tolerance may not be achievable if the quadrature rule is not accurate enough, or if τ is not sufficiently small, and therefore a maximum number of iterations maxiter is provided to avoid an infinite loops. In such cases, it is recommended that you either increase the quadrature order, or decrease the mesh size.

Note: this is obviously a heuristic, and may not be accurate in all cases.

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.

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.

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)
meshgen(Ω, n_dict)
meshgen(Ω; meshsize)

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.

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

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.

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 derivative of the single- and double-layer, respectively).

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.
• (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 thesourcecurve/surface (which is assumed to be closed). Whentarget === source,target_location` is not needed.

single_double_layer_potential(; op, source)

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

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 [8] 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.