Add integral and localintegral

Project.toml

@@ -11,6 +11,7 @@ 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"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
@@ -38,6 +39,7 @@ CoordRefSystems = "0.19"
 DelaunayTriangulation = "1.0"
 DifferentiationInterface = "0.7"
 Distances = "0.10"
+FastGaussQuadrature = "1.1"
 FiniteDifferences = "0.12"
 LinearAlgebra = "1.10"
 Makie = "0.24"

src/Meshes.jl

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

src/integration.jl

@@ -0,0 +1,45 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    integral(fun, geom; order=1)
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+[`Point`](@ref)s to values in a linear space.
+
+The quadrature `order` can be specified to increase the accuracy.
+
+See also [`localintegral`](@ref).
+"""
+integral(fun, geom; order=1) = localintegral(fun ∘ geom, geom; order)
+
+"""
+    localintegral(fun, geom; order=1)
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+parametric coordinates `uvw` to values in a linear space.
+
+The quadrature `order` can be specified to increase the accuracy.
+
+See also [`integral`](@ref).
+"""
+function localintegral(fun, geom; order=1)
+  # parametric dimension
+  dim = paramdim(geom)
+
+  # Gauss-Legendre quadrature points and weights
+  ts, ws = gausslegendre(order)
+  tgrid = Iterators.product(ntuple(_ -> ts, dim)...)
+  wgrid = Iterators.product(ntuple(_ -> ws, dim)...)
+
+  # compute integral with change of variable and differential element
+  Σwᵢfᵢ = sum(zip(tgrid, wgrid)) do (t, w)
+    # change of variable [-1, 1] → [0, 1]
+    uvw = ntuple(i -> (t[i] + 1) / 2, dim)
+    prod(w) * fun(uvw...) * differential(geom, uvw)
+  end
+
+  # adjust for change of variable
+  Σwᵢfᵢ / 2^dim
+end

test/integration.jl

@@ -0,0 +1,5 @@
+@testitem "integral" setup = [Setup] begin
+end
+
+@testitem "localintegral" setup = [Setup] begin
+end