Add support for integrating over PolyArea and Domain types

CHANGELOG.md

@@ -7,6 +7,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- Adds integral methods for `PolyArea` and subtypes of `Domain` with explicit tests for `CartesianGrid`, `PolyArea`, `RegularGrid`, `SimpleMesh`, and `StructuredGrid`.
+
+### Changed
+
+- Generalizes the alias functions (e.g. `lineintegral`) to also accept `Domain`s.
+
+
 ## [0.16.3] - 2025-06-26
 
 ### Changed

docs/src/support.md

@@ -26,14 +26,15 @@ The following Support Matrix captures the current state of support for all geome
 combinations. Entries with a green check mark are fully supported and pass unit tests
 designed to check for accuracy.
 
-| `Meshes.Geometry` | `GaussKronrod` | `GaussLegendre` | `HAdaptiveCubature` |
+| `Meshes.Geometry/Domain` | `GaussKronrod` | `GaussLegendre` | `HAdaptiveCubature` |
 |----------|----------------|---------------|---------------------|
 | `Ball` in `𝔼{2}` | ⚠️ | ✅ | ✅ |
 | `Ball` in `𝔼{3}` | 🛑 | ✅ | ✅ |
 | `BezierCurve` | ✅ | ✅ | ✅ |
 | `Box` in `𝔼{1}` | ✅ | ✅ | ✅ |
 | `Box` in `𝔼{2}` | ⚠️ | ✅ | ✅ |
 | `Box` in `𝔼{≥3}` | 🛑 | ✅ | ✅ |
+| `CartesianGrid` | ✅ | ✅ | ✅ |
 | `Circle` | ✅ | ✅ | ✅ |
 | `Cone` | 🛑 | ✅ | ✅ |
 | `ConeSurface` | ⚠️ | ✅ | ✅ |
@@ -48,16 +49,18 @@ designed to check for accuracy.
 | `ParaboloidSurface` | ⚠️ | ✅ | ✅ |
 | `ParametrizedCurve` | ✅ | ✅ | ✅ |
 | `Plane` | ✅ | ✅ | ✅ |
-| `Polyarea` | 🛑 | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
+| `PolyArea` | ⚠️ | ✅ | ✅ |
 | `Pyramid` | ⚠️ | ✅ | ✅ |
 | `Quadrangle` | ⚠️ | ✅ | ✅ |
 | `Ray` | ✅ | ✅ | ✅ |
+| `RegularGrid` | ✅ | ✅ | ✅ |
 | `Ring` | ✅ | ✅ | ✅ |
 | `Rope` | ✅ | ✅ | ✅ |
 | `Segment` | ✅ | ✅ | ✅ |
-| `SimpleMesh` | 🛑 | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) |
+| `SimpleMesh` | ⚠️ | ✅ | ✅ |
 | `Sphere` in `𝔼{2}` | ✅ | ✅ | ✅ |
 | `Sphere` in `𝔼{3}` | ⚠️ | ✅ | ✅ |
+| `StructuredGrid` | ✅ | ✅ | ✅ |
 | `Tetrahedron` | ⚠️ | ✅ | ✅ |
 | `Triangle` | ✅ | ✅ | ✅ |
 | `Torus` | ⚠️ | ✅ | ✅ |
@@ -66,6 +69,5 @@ designed to check for accuracy.
 | Symbol | Support Level |
 |--------|---------|
 | ✅ | Supported |
-| 🎗️ | Planned to support in the future |
 | ⚠️ | Deprecated |
-| 🛑 | Not supported |
+| 🛑 | Not supported |

ext/MeshIntegralsEnzymeExt.jl

@@ -1,13 +1,13 @@
 module MeshIntegralsEnzymeExt
 
-using MeshIntegrals: MeshIntegrals, AutoEnzyme
-using Meshes: Meshes
-using Enzyme: Enzyme
+import MeshIntegrals
+import Meshes
+import Enzyme
 
 function MeshIntegrals.jacobian(
         geometry::Meshes.Geometry,
         ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
-        ::AutoEnzyme
+        ::MeshIntegrals.AutoEnzyme
 ) where {T <: AbstractFloat}
     Dim = Meshes.paramdim(geometry)
     if Dim != length(ts)
@@ -16,9 +16,10 @@ function MeshIntegrals.jacobian(
     return Meshes.to.(Enzyme.jacobian(Enzyme.Forward, geometry, ts...))
 end
 
-# Supports all geometries except for those that throw errors
-# See GitHub Issue #154 for more information
+# Supports all geometries/domains by default
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.Geometry}) = true
+MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.Domain}) = true
+# except for those that throw errors (see GitHub Issue #154)
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.BezierCurve}) = false
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.CylinderSurface}) = false
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.Cylinder}) = false

src/MeshIntegrals.jl

@@ -9,6 +9,9 @@ import LinearAlgebra
 import QuadGK
 import Unitful
 
+# Same as the internal-only Meshes.GeometryOrDomain
+const GeometryOrDomain = Union{Meshes.Geometry, Meshes.Domain}
+
 include("differentiation.jl")
 export DifferentiationMethod, FiniteDifference, AutoEnzyme, jacobian
 
@@ -31,6 +34,7 @@ include("specializations/CylinderSurface.jl")
 include("specializations/FrustumSurface.jl")
 include("specializations/Line.jl")
 include("specializations/Plane.jl")
+include("specializations/PolyArea.jl")
 include("specializations/Ray.jl")
 include("specializations/Ring.jl")
 include("specializations/Rope.jl")

src/integral.jl

@@ -32,6 +32,18 @@ function integral(
     _integral(f, geometry, rule; kwargs...)
 end
 
+function integral(
+        f,
+        domain::Meshes.Domain,
+        rule::I = Meshes.paramdim(domain) == 1 ? GaussKronrod() : HAdaptiveCubature();
+        kwargs...
+) where {I <: IntegrationRule}
+    # Discretize the Domain into primitive geometries, sum the integrals over those
+    subintegral(geometry) = integral(f, geometry, rule; kwargs...)
+    subgeometries = Meshes.elements(Meshes.discretize(domain)) |> collect
+    return sum(subintegral, subgeometries)
+end
+
 ################################################################################
 #                            Integral Workers
 ################################################################################

src/integral_aliases.jl

@@ -16,7 +16,7 @@ Rule types available:
 """
 function lineintegral(
         f,
-        geometry::Geometry,
+        geometry::GeometryOrDomain,
         rule::IntegrationRule = GaussKronrod();
         kwargs...
 )
@@ -48,7 +48,7 @@ Algorithm types available:
 """
 function surfaceintegral(
         f,
-        geometry::Geometry,
+        geometry::GeometryOrDomain,
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
 )
@@ -80,7 +80,7 @@ Algorithm types available:
 """
 function volumeintegral(
         f,
-        geometry::Geometry,
+        geometry::GeometryOrDomain,
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
 )

src/specializations/PolyArea.jl

@@ -0,0 +1,60 @@
+################################################################################
+#                      Specialized Methods for PolyArea
+#
+# Why Specialized?
+#   The PolyArea geometry defines a surface bounded by a polygon, but Meshes.jl
+#   cannot define a parametric function for a polyarea because a general solution
+#   for one does not exist. However, they can be partitioned into simpler elements
+#   and integrated separately before finally summing the result.
+################################################################################
+
+"""
+    integral(f, area::PolyArea, rule = HAdaptiveCubature();
+             diff_method=FiniteDifference(), FP=Float64)
+
+Like [`integral`](@ref) but integrates over the surface domain defined by a `PolyArea`.
+The surface is first discretized into facets that are integrated independently using
+the specified integration `rule`.
+
+# Arguments
+- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
+- `area`: a `PolyArea` that defines the integration domain
+- `rule = HAdaptiveCubature()`: optionally, the `IntegrationRule` used for integration
+
+# Keyword Arguments
+- `diff_method::DifferentiationMethod`: the method to use for
+calculating Jacobians that are used to calculate differential elements
+- `FP = Float64`: the floating point precision desired
+"""
+function integral(
+        f,
+        area::Meshes.PolyArea,
+        rule::I;
+        kwargs...
+) where {I <: IntegrationRule}
+    # Partition the PolyArea, sum the integrals over each of those areas
+    subintegral(area) = integral(f, area, rule; kwargs...)
+    subgeometries = collect(Meshes.elements(Meshes.discretize(area)))
+    return sum(subintegral, subgeometries)
+end
+
+"""
+    surfaceintegral(f, area, rule = HAdaptiveCubature(); FP = Float64)
+
+Numerically integrate a given function `f(::Point)` over the surface domain defined
+by a `PolyArea`. The surface is first discretized into facets that are integrated
+independently using the specified integration `rule`.
+
+Algorithm types available:
+- [`GaussKronrod`](@ref)
+- [`GaussLegendre`](@ref)
+- [`HAdaptiveCubature`](@ref) (default)
+"""
+function surfaceintegral(
+        f,
+        area::Meshes.PolyArea,
+        rule::IntegrationRule = HAdaptiveCubature();
+        kwargs...
+)
+    return integral(f, area, rule; kwargs...)
+end

src/specializations/Ring.jl

@@ -33,5 +33,6 @@ function integral(
         kwargs...
 ) where {I <: IntegrationRule}
     # Convert the Ring into Segments, sum the integrals of those 
-    return sum(segment -> _integral(f, segment, rule; kwargs...), Meshes.segments(ring))
+    subintegral(segment) = _integral(f, segment, rule; kwargs...)
+    return sum(subintegral, Meshes.segments(ring))
 end

src/specializations/Rope.jl

@@ -33,5 +33,6 @@ function integral(
         kwargs...
 ) where {I <: IntegrationRule}
     # Convert the Rope into Segments, sum the integrals of those
-    return sum(segment -> integral(f, segment, rule; kwargs...), Meshes.segments(rope))
+    subintegral(segment) = integral(f, segment, rule; kwargs...)
+    return sum(subintegral, Meshes.segments(rope))
 end

src/utils.jl

@@ -14,7 +14,7 @@ end
 
 """
     supports_autoenzyme(geometry::Geometry)
-    supports_autoenzyme(type::Type{<:Geometry})
+    supports_autoenzyme(type::Type{<:GeometryOrDomain})
 
 Return whether a geometry or geometry type has a parametric function that can be
 differentiated with Enzyme. See GitHub Issue #154 for more information.
@@ -25,7 +25,7 @@ function supports_autoenzyme end
 supports_autoenzyme(::Type{<:Any}) = false
 
 # If provided a geometry instance, re-run with the type as argument
-supports_autoenzyme(::G) where {G <: Geometry} = supports_autoenzyme(G)
+supports_autoenzyme(::G) where {G <: GeometryOrDomain} = supports_autoenzyme(G)
 
 """
     _check_diff_method_support(::Geometry, ::DifferentiationMethod) -> nothing