API

guiggiani_singular_integral(
	K::Inti.AbstractKernel,
	û,
	x̂,
	el::Inti.ReferenceInterpolant,
	quad_rho,
	quad_theta,
	method::AbstractMethod = FullRichardsonExpansion(),
)

Compute the singular integral of kernel K over element el using Guiggiani's method.

This function evaluates:

∫_{el} K(x, y) * û(ŷ) dS(y)

where x = el(x̂) is the singular point on the element.

The method works by:

1. Transforming to polar coordinates centered at x̂
2. Computing Laurent coefficients (f₋₂, f₋₁) to extract the singularity
3. Integrating the regularized integrand: K - f₋₂/ρ² - f₋₁/ρ
4. Adding back the analytical integrals of the singular terms

Arguments

• K::Inti.AbstractKernel: The kernel (must have a defined singularity_order)
• û: Density function defined on the reference element
• x̂: Singular point location on the reference element
• el::Inti.ReferenceInterpolant: The reference element
• quad_rho: Quadrature rule for the radial direction (e.g., Inti.GaussLegendre(10))
• quad_theta: Quadrature rule for the angular direction (e.g., Inti.GaussLegendre(20))
• method::AbstractMethod: Method for computing Laurent coefficients (see laurents_coeffs)

Returns

• The value of the singular integral

Notes

• The singularity order is automatically detected from Inti.singularity_order(K)
• In polar coordinates, the order is adjusted by +1 (due to the Jacobian factor ρ)
• Currently supports singularity orders -1, -2, and -3 in Cartesian coordinates

Examples

# Basic usage with default parameters
K = Inti.Laplace(dim=2) |> Inti.HyperSingularKernel()
el = # ... reference element ...
û = ŷ -> 1.0  # constant density
x̂ = SVector(0.5, 0.5)
quad_rho = Inti.GaussLegendre(10)
quad_theta = Inti.GaussLegendre(20)

I = guiggiani_singular_integral(K, û, x̂, el, quad_rho, quad_theta)

# Using automatic differentiation for speed
I = guiggiani_singular_integral(K, û, x̂, el, quad_rho, quad_theta, AutoDiffExpansion())

# With custom Richardson parameters
params = RichardsonParams(atol=1e-12, maxeval=10)
I = guiggiani_singular_integral(K, û, x̂, el, quad_rho, quad_theta, FullRichardsonExpansion(params))

laurents_coeffs(
	K::Inti.AbstractKernel, 
	el::Inti.ReferenceInterpolant, 
	û, 
	x̂,
	method::AbstractMethod = FullRichardsonExpansion(),
)

Compute Laurent coefficients (f₋₂, f₋₁) for the kernel K in polar coordinates centered at x̂.

Arguments

• K::Inti.AbstractKernel: The kernel to expand
• el::Inti.ReferenceInterpolant: The reference element
• û: Function defined on the reference element (density function)
• x̂: Point on the reference element (singularity location)
• method::AbstractMethod: Expansion method to use. Can be:
  • AnalyticalExpansion(): Uses closed-form analytical expressions (fastest, limited availability)
  • AutoDiffExpansion(): Uses automatic differentiation (fast, requires translation-invariant kernels)
  • SemiRichardsonExpansion(params): Hybrid Richardson extrapolation (moderate speed)
  • FullRichardsonExpansion(params): Full Richardson extrapolation (slowest, always available)

Returns

• ℒ(θ): A memoized function that returns (f₋₂, f₋₁) for a given angle θ

Examples

# Using default method (FullRichardson)
ℒ = laurents_coeffs(K, el, ori, û, x̂)
f₋₂, f₋₁ = ℒ(0.5)  # Evaluate at θ = 0.5

# Using AutoDiff method
ℒ = laurents_coeffs(K, el, ori, û, x̂, AutoDiffExpansion())

# Using explicit method with custom parameters
params = RichardsonParams(atol=1e-10, rtol=1e-8, maxeval=10)
ℒ = laurents_coeffs(K, el, ori, û, x̂, FullRichardsonExpansion(params))