Add integrals over geometries

Project.toml

@@ -11,7 +11,9 @@ CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+IntegrationInterface = "03ca86a4-b7cd-4b10-a38f-d7a84f40de1e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -38,7 +40,9 @@ CoordRefSystems = "0.19"
 DelaunayTriangulation = "1.0"
 DifferentiationInterface = "0.7"
 Distances = "0.10"
+FastGaussQuadrature = "1.1"
 FiniteDifferences = "0.12"
+IntegrationInterface = "0.1"
 LinearAlgebra = "1.10"
 Makie = "0.24"
 NearestNeighbors = "0.4"

src/Meshes.jl

@@ -44,6 +44,10 @@ import NearestNeighbors: MinkowskiMetric
 import DifferentiationInterface as DI
 import FiniteDifferences: central_fdm
 
+# Integration API
+import IntegrationInterface as II
+import FastGaussQuadrature: gausslegendre
+
 # Transforms API
 import TransformsBase: Transform, →
 import TransformsBase: isrevertible, isinvertible
@@ -104,6 +108,7 @@ include("predicates.jl")
 
 # calculus
 include("differentation.jl")
+include("integration.jl")
 
 # operations
 include("centroid.jl")
@@ -458,6 +463,8 @@ export
   derivative,
   jacobian,
   differential,
+  integral,
+  localintegral,
 
   # centroids
   centroid,

src/integration.jl

@@ -0,0 +1,97 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# default integration method
+const GAUSSLEGENDRE = II.Backend.Quadrature(gausslegendre(20))
+
+"""
+    integral(fun, geom[, method]) 
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+[`Point`](@ref)s to values in a linear space.
+
+    integral(fun, dom[, method])
+
+Alternatively, calculate the integral over the `dom`ain (e.g., mesh) by
+summing the integrals for each constituent geometry.
+
+By default, use Gauss-Legendre quadrature rule with `n` nodes so that
+polynomials of degree up to `2n-1` are integrated exactly.
+
+See also [`localintegral`](@ref).
+"""
+integral(fun, geom::Geometry, method=GAUSSLEGENDRE) = localintegral(fun ∘ geom, geom, method)
+
+# cylinder surface is the union of the curved surface and the top and bottom disks
+integral(fun, cylsurf::CylinderSurface, method=GAUSSLEGENDRE) =
+  localintegral(fun ∘ cylsurf, cylsurf, method) +
+  integral(fun, top(cylsurf), method) +
+  integral(fun, bottom(cylsurf), method)
+
+# cone surface is the union of the curved surface and the base disk
+integral(fun, conesurf::ConeSurface, method=GAUSSLEGENDRE) =
+  localintegral(fun ∘ conesurf, conesurf, method) + integral(fun, base(conesurf), method)
+
+# frustum surface is the union of the curved surface and the top and bottom disks
+integral(fun, frustumsurf::FrustumSurface, method=GAUSSLEGENDRE) =
+  localintegral(fun ∘ frustumsurf, frustumsurf, method) +
+  integral(fun, top(frustumsurf), method) +
+  integral(fun, bottom(frustumsurf), method)
+
+# multi-geometry is the union of its constituent geometries
+integral(fun, multi::Multi, method=GAUSSLEGENDRE) = sum(integral(fun, geom, method) for geom in parent(multi))
+
+# domain is the union of its constituent geometries
+integral(fun, dom::Domain, method=GAUSSLEGENDRE) = sum(integral(fun, geom, method) for geom in dom)
+
+"""
+    localintegral(fun, geom[, method])
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+parametric coordinates `uvw` to values in a linear space.
+
+By default, use Gauss-Legendre quadrature rule with `n` nodes so that
+polynomials of degree up to `2n-1` are integrated exactly.
+
+See also [`integral`](@ref).
+"""
+function localintegral(fun, geom::Geometry, method=GAUSSLEGENDRE)
+  # integrand is equal to function times differential element
+  integrand(uvw) = fun(uvw...) * differential(geom, uvw)
+
+  # domain of integration for the given geometry
+  domain = ∫domain(geom)
+
+  # perform numerical integration
+  II.integral(integrand, domain; backend=method)()
+end
+
+function ∫domain(geom::Geometry)
+  N = paramdim(geom)
+  T = numtype(lentype(geom))
+  a = ntuple(_ -> zero(T), N)
+  b = ntuple(_ -> one(T), N)
+  II.Domain.Box(a, b)
+end
+
+function ∫domain(ray::Ray)
+  T = numtype(lentype(ray))
+  a = (zero(T),)
+  b = (II.Infinity(one(T)),)
+  II.Domain.Box(a, b)
+end
+
+function ∫domain(line::Line)
+  T = numtype(lentype(line))
+  a = (-II.Infinity(one(T)),)
+  b = (II.Infinity(one(T)),)
+  II.Domain.Box(a, b)
+end
+
+function ∫domain(plane::Plane)
+  T = numtype(lentype(plane))
+  a = (-II.Infinity(one(T)), -II.Infinity(one(T)))
+  b = (II.Infinity(one(T)), II.Infinity(one(T)))
+  II.Domain.Box(a, b)
+end