Add support for integrating over PolyArea and Domain types (#182)

* Add tests for SimpleMesh

* Add integral method for Domains, generalize aliases to GeometryOrDomain

* Define specialization methods for PolyArea

* Use explicit namespace for Meshes.Domain

* Generalize support_status geometry type

* Bugfix, explicit import GeometryOrDomain

* Explicit use of non-exported union

* Explicit namespace bugfix

* Extend supports_autoenzyme to Domains

* Formatting changes

* Adjust test imports

* Add tests for RegularGrid

* Add tests for CartesianGrid, add missing ustrip import

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Add explicit rtol's, reduce dims to speed up RegularGrid tests

* Generalize supports_autoenzyme from utils to Domains

* Bugfix update old name

* Make partitioned integral codes consistent

* Add PolyArea tests

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Bugfix

* Bugfix

* Reverse hole region coord path

* Apply suggestion

Co-authored-by: Joshua Lampert <[email protected]>

* Mirror GeometryOrDomain internally

* Bugfix

* Remove unnecessary use of GeometryOrDomain, switch usings to imports

* Update docs and merge recent v0.16.3 release changes

* Add tests and docs for StructuredGrid

* Update test/combinations.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Remove unnecessary method

* Docstring and style updates

* Update docstrings for consistency and clarity, associated style improvements

* Add note to changelog

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <[email protected]>

* Switch to function format due to line length

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Joshua Lampert <[email protected]>

Author: 3 people

CHANGELOG.md

@@ -7,6 +7,16 @@ 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.
+- Updated docstrings to improve clarity and consistency.
+
+
 ## [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

@@ -3,35 +3,52 @@
 ################################################################################
 
 """
-    integral(f, geometry[, rule]; diff_method=_default_diff_method(geometry, FP), FP=Float64)
+    integral(f, geometry[, rule]; kwargs...)
 
 Numerically integrate a given function `f(::Point)` over the domain defined by
 a `geometry` using a particular numerical integration `rule` with floating point
 precision of type `FP`.
 
 # Arguments
 - `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
-- `geometry`: some `Meshes.Geometry` that defines the integration domain
+- `geometry`: a Meshes.jl `Geometry` or `Domain` that defines the integration domain
 - `rule`: optionally, the `IntegrationRule` used for integration (by default
 `GaussKronrod()` in 1D and `HAdaptiveCubature()` else)
 
 # Keyword Arguments
-- `diff_method::DifferentiationMethod = _default_diff_method(geometry, FP)`: the method to
-use for calculating Jacobians that are used to calculate differential elements
-- `FP = Float64`: the floating point precision desired.
+- `diff_method::DifferentiationMethod`: manually specifies the differentiation method use to
+calculate Jacobians within the integration domain.
+- `FP = Float64`: manually specifies the desired output floating point precision
 """
 function integral end
 
+# Default integration rule to use if unspecified
+function _default_rule(geometry)
+    ifelse(Meshes.paramdim(geometry) == 1, GaussKronrod(), HAdaptiveCubature())
+end
+
 # If only f and geometry are specified, select default rule
 function integral(
         f,
         geometry::Geometry,
-        rule::I = Meshes.paramdim(geometry) == 1 ? GaussKronrod() : HAdaptiveCubature();
+        rule::I = _default_rule(geometry);
         kwargs...
 ) where {I <: IntegrationRule}
     _integral(f, geometry, rule; kwargs...)
 end
 
+function integral(
+        f,
+        domain::Meshes.Domain,
+        rule::I = _default_rule(domain);
+        kwargs...
+) where {I <: IntegrationRule}
+    # Discretize the Domain into primitive geometries, sum the integrals over those
+    subintegral(geometry) = integral(f, geometry, rule; kwargs...)
+    subgeometries = collect(Meshes.elements(Meshes.discretize(domain)))
+    return sum(subintegral, subgeometries)
+end
+
 ################################################################################
 #                            Integral Workers
 ################################################################################

src/integral_aliases.jl

@@ -9,25 +9,26 @@ Numerically integrate a given function `f(::Point)` along a line-like `geometry`
 using a particular numerical integration `rule` with floating point precision of
 type `FP`.
 
+This is a convenience wrapper around [`integral`](@ref) that additionally enforces
+a requirement that the geometry have one parametric dimension.
+
 Rule types available:
 - [`GaussKronrod`](@ref) (default)
 - [`GaussLegendre`](@ref)
 - [`HAdaptiveCubature`](@ref)
 """
 function lineintegral(
         f,
-        geometry::Geometry,
+        geometry::GeometryOrDomain,
         rule::IntegrationRule = GaussKronrod();
         kwargs...
 )
-    N = Meshes.paramdim(geometry)
-
-    if N == 1
-        return integral(f, geometry, rule; kwargs...)
-    else
+    if (N = Meshes.paramdim(geometry)) != 1
         throw(ArgumentError("Performing a line integral on a geometry \
                             with $N parametric dimensions not supported."))
     end
+
+    return integral(f, geometry, rule; kwargs...)
 end
 
 ################################################################################
@@ -41,25 +42,26 @@ Numerically integrate a given function `f(::Point)` along a surface `geometry`
 using a particular numerical integration `rule` with floating point precision of
 type `FP`.
 
+This is a convenience wrapper around [`integral`](@ref) that additionally enforces
+a requirement that the geometry have two parametric dimensions.
+
 Algorithm types available:
 - [`GaussKronrod`](@ref)
 - [`GaussLegendre`](@ref)
 - [`HAdaptiveCubature`](@ref) (default)
 """
 function surfaceintegral(
         f,
-        geometry::Geometry,
+        geometry::GeometryOrDomain,
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
 )
-    N = Meshes.paramdim(geometry)
-
-    if N == 2
-        return integral(f, geometry, rule; kwargs...)
-    else
+    if (N = Meshes.paramdim(geometry)) != 2
         throw(ArgumentError("Performing a surface integral on a geometry \
                             with $N parametric dimensions not supported."))
     end
+
+    return integral(f, geometry, rule; kwargs...)
 end
 
 ################################################################################
@@ -73,23 +75,24 @@ Numerically integrate a given function `f(::Point)` throughout a volumetric
 `geometry` using a particular numerical integration `rule` with floating point
 precision of type `FP`.
 
+This is a convenience wrapper around [`integral`](@ref) that additionally enforces
+a requirement that the geometry have three parametric dimensions.
+
 Algorithm types available:
 - [`GaussKronrod`](@ref)
 - [`GaussLegendre`](@ref)
 - [`HAdaptiveCubature`](@ref) (default)
 """
 function volumeintegral(
         f,
-        geometry::Geometry,
+        geometry::GeometryOrDomain,
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
 )
-    N = Meshes.paramdim(geometry)
-
-    if N == 3
-        return integral(f, geometry, rule; kwargs...)
-    else
+    if (N = Meshes.paramdim(geometry)) != 3
         throw(ArgumentError("Performing a volume integral on a geometry \
                             with $N parametric dimensions not supported."))
     end
+
+    return integral(f, geometry, rule; kwargs...)
 end

src/specializations/BezierCurve.jl

@@ -13,20 +13,11 @@
 #                              integral
 ################################################################################
 """
-    integral(f, curve::BezierCurve, rule = GaussKronrod();
-             diff_method=Analytical(), FP=Float64, alg=Meshes.Horner())
+    integral(f, curve::BezierCurve[, rule = GaussKronrod()]; kwargs...)
 
 Like [`integral`](@ref) but integrates along the domain defined by `curve`.
 
-# Arguments
-- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
-- `curve`: a `Meshes.BezierCurve` that defines the integration domain
-- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
-
-# Keyword Arguments
-- `diff_method::DifferentiationMethod = Analytical()`: the method to use for
-calculating Jacobians that are used to calculate differential elements
-- `FP = Float64`: the floating point precision desired
+# Special Keyword Arguments
 - `alg = Meshes.Horner()`:  the method to use for parametrizing `curve`. Alternatively,
 `alg=Meshes.DeCasteljau()` can be specified for increased accuracy, but comes with a
 steep performance cost, especially for curves with a large number of control points.

src/specializations/PolyArea.jl

@@ -0,0 +1,28 @@
+################################################################################
+#                      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()]; kwargs...)
+
+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`.
+"""
+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

src/specializations/Ring.jl

@@ -9,22 +9,11 @@
 ################################################################################
 
 """
-    integral(f, ring::Ring, rule = GaussKronrod();
-             diff_method=FiniteDifference(), FP=Float64)
+    integral(f, ring::Ring[, rule = GaussKronrod()]; kwargs...)
 
 Like [`integral`](@ref) but integrates along the domain defined by `ring`. The
 specified integration `rule` is applied independently to each segment formed by
 consecutive points in the Ring.
-
-# Arguments
-- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
-- `ring`: a `Ring` that defines the integration domain
-- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
-
-# Keyword Arguments
-- `diff_method::DifferentiationMethod = FiniteDifference()`: the method to use for
-calculating Jacobians that are used to calculate differential elements
-- `FP = Float64`: the floating point precision desired
 """
 function integral(
         f,
@@ -33,5 +22,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

@@ -9,22 +9,11 @@
 ################################################################################
 
 """
-    integral(f, rope::Rope, rule = GaussKronrod();
-             diff_method=FiniteDifference(), FP=Float64)
+    integral(f, rope::Rope[, rule = GaussKronrod()]; kwargs...)
 
 Like [`integral`](@ref) but integrates along the domain defined by `rope`. The
 specified integration `rule` is applied independently to each segment formed by
 consecutive points in the Rope.
-
-# Arguments
-- `f`: an integrand function, i.e. any callable with a method `f(::Meshes.Point)`
-- `rope`: a `Rope` that defines the integration domain
-- `rule = GaussKronrod()`: optionally, the `IntegrationRule` used for integration
-
-# Keyword Arguments
-- `diff_method::DifferentiationMethod = FiniteDifference()`: the method to use for
-calculating Jacobians that are used to calculate differential elements
-- `FP = Float64`: the floating point precision desired
 """
 function integral(
         f,
@@ -33,5 +22,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