MeshIntegrals.jl uses differential forms to enable fast and easy numerical integration of arbitrary integrand functions over domains defined via Meshes.jl geometries. This is achieved using:

• Gauss-Legendre quadrature rules from FastGaussQuadrature.jl: GaussLegendre(n)
• H-adaptive Gauss-Kronrod quadrature rules from QuadGK.jl: GaussKronrod(kwargs...)
• H-adaptive cubature rules from HCubature.jl: HAdaptiveCubature(kwargs...)

These solvers have support for integrand functions that produce scalars, vectors, and Unitful.jl Quantity types. While HCubature.jl does not natively support Quantity type integrands, this package provides a compatibility layer to enable this feature.

integral(f, geometry)

Performs a numerical integration of some integrand function f(p::Meshes.Point) over the domain specified by geometry. A default integration method will be automatically selected according to the geometry: GaussKronrod() for 1D, and HAdaptiveCubature() for all others.

integral(f, geometry, algorithm, FP=Float64)

Performs a numerical integration of some integrand function f(p::Meshes.Point) over the domain specified by geometry using the specified integration algorithm, e.g. GaussKronrod().

Optionally, a fourth argument can be provided to specify the floating point precision level desired. This setting can be manipulated if your integrand function produces outputs with alternate floating point precision (e.g. Float16, BigFloat, etc) AND you'd prefer to avoid implicit type promotions.

lineintegral(f, geometry)
surfaceintegral(f, geometry)
volumeintegral(f, geometry)

Alias functions are provided for convenience. These are simply wrappers for integral that first validate that the provided geometry has the expected number of parametric/manifold dimensions. Like with integral in the examples above, the algorithm can also be specified as a third-argument.

• lineintegral for curve-like geometries or polytopes (e.g. Segment, Ray, BezierCurve, Rope, etc)
• surfaceintegral for surfaces (e.g. Disk, Sphere, CylinderSurface, etc)
• volumeintegral for (3D) volumes (e.g. Ball, Cone, Torus, etc)

using Meshes
using MeshIntegrals
using Unitful

# Define a path that approximates a sine-wave on the xy-plane
mypath = BezierCurve(
    [Point(t*u"m", sin(t)*u"m", 0.0u"m") for t in range(-pi, pi, length=361)]
)

# Map f(::Point) -> f(x, y, z) in unitless coordinates
f(p::Meshes.Point) = f(ustrip(to(p))...)

# Integrand function in units of Ohms/meter
f(x, y, z) = (1 / sqrt(1 + cos(x)^2)) * u"Ω/m"

integral(f, mypath)
# -> Approximately 2*Pi Ω