Add support for integrating Meshes.Domain and PolyArea

ext/MeshIntegralsEnzymeExt.jl

@@ -18,7 +18,7 @@ end
 
 # Supports all geometries except for those that throw errors
 # See GitHub Issue #154 for more information
-MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.Geometry}) = true
+MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.GeometryOrDomain}) = true
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.BezierCurve}) = false
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.CylinderSurface}) = false
 MeshIntegrals.supports_autoenzyme(::Type{<:Meshes.Cylinder}) = false

src/MeshIntegrals.jl

@@ -1,7 +1,7 @@
 module MeshIntegrals
 using CliffordNumbers: CliffordNumbers, VGA, ∧
 using CoordRefSystems: CoordRefSystems, CRS
-using Meshes: Meshes, Geometry
+using Meshes: Meshes, Geometry, GeometryOrDomain
 
 import FastGaussQuadrature
 import HCubature
@@ -31,6 +31,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))
+    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,59 @@
+################################################################################
+#                      Specialized Methods for PolyArea
+#
+# Why Specialized?
+#   The PolyArea geometry defines a surface bounded by a polygon but Meshes.jl
+#   can not 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...)
+    return sum(subintegral, Meshes.discretize(area))
+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 <: Meshes.GeometryOrDomain} = supports_autoenzyme(G)
 
 """
     _check_diff_method_support(::Geometry, ::DifferentiationMethod) -> nothing

test/combinations.jl

@@ -15,10 +15,10 @@ This file includes tests for:
 
 @testsnippet Combinations begin
     using CoordRefSystems
-    using LinearAlgebra: norm
+    import LinearAlgebra: norm
     using Meshes
     using MeshIntegrals
-    using Unitful
+    import Unitful: @u_str, Quantity, ustrip
     import Enzyme
 
     # Used for testing callable objects as integrand functions
@@ -28,7 +28,7 @@ This file includes tests for:
     (c::Callable)(p::Meshes.Point) = c.f(p)
 
     # Stores a testable combination
-    struct TestableGeometry{F <: Function, G <: Geometry, U <: Unitful.Quantity}
+    struct TestableGeometry{F <: Function, G <: Meshes.GeometryOrDomain, U <: Quantity}
         integrand::F
         geometry::G
         solution::U
@@ -49,7 +49,7 @@ This file includes tests for:
     end
 
     # Shortcut constructor for geometries with typical support structure
-    function SupportStatus(geometry::G;) where {G <: Geometry}
+    function SupportStatus(geometry::G) where {G <: Meshes.GeometryOrDomain}
         # Check whether AutoEnzyme should be supported, i.e. not on blacklist
         unsupported_Gs = Union{BezierCurve, Cylinder, CylinderSurface, ParametrizedCurve}
         autoenzyme = !(G <: unsupported_Gs)
@@ -272,6 +272,26 @@ end
     runtests(testable; rtol = 1e-6)
 end
 
+@testitem "Meshes.CartesianGrid" setup=[Combinations] begin
+    # Geometry
+    a = π
+    start = Point(0, 0)
+    finish = Point(a, a)
+    dims = (4, 4)
+    grid = CartesianGrid(start, finish, dims = dims)
+
+    # Integrand & Solution
+    function integrand(p::Meshes.Point)
+        x₁, x₂ = ustrip.(to(p))
+        (√(a^2 - x₁^2) + √(a^2 - x₂^2)) * u"A"
+    end
+    solution = 2a * (π * a^2 / 4) * u"A*m^2"
+
+    # Package and run tests
+    testable = TestableGeometry(integrand, grid, solution)
+    runtests(testable; rtol = 1e-6)
+end
+
 @testitem "Meshes.Circle" setup=[Combinations] begin
     # Geometry
     center = Point(1, 2, 3)
@@ -628,6 +648,26 @@ end
     runtests(testable)
 end
 
+@testitem "Meshes.RegularGrid" setup=[Combinations] begin
+    # Geometry
+    a = π
+    start = Point(0, 0)
+    finish = Point(a, a)
+    dims = (4, 4)
+    grid = RegularGrid(start, finish, dims = dims)
+
+    # Integrand & Solution
+    function integrand(p::Meshes.Point)
+        x₁, x₂ = ustrip.(to(p))
+        (√(a^2 - x₁^2) + √(a^2 - x₂^2)) * u"A"
+    end
+    solution = 2a * (π * a^2 / 4) * u"A*m^2"
+
+    # Package and run tests
+    testable = TestableGeometry(integrand, grid, solution)
+    runtests(testable; rtol = 1e-6)
+end
+
 @testitem "Meshes.Ring" setup=[Combinations] begin
     # Geometry
     a = Point(0, 0, 0)
@@ -677,7 +717,7 @@ end
 
     # Integrand & Solution
     a, b = (7.1, 4.6)  # arbitrary constants > 0
-    function integrand(p::P; a = a, b = b) where {P <: Meshes.Point}
+    function integrand(p::Meshes.Point; a = a, b = b)
         r = ustrip(u"m", norm(to(p)))
         exp(r * log(a) + (1 - r) * log(b)) * u"A"
     end
@@ -688,6 +728,26 @@ end
     runtests(testable)
 end
 
+@testitem "Meshes.SimpleMesh" setup=[Combinations] begin
+    # Geometry
+    a = π
+    points = [(0, 0), (a, 0), (0, a), (a, a), (0.25a, 0.5a), (0.75a, 0.5a)]
+    tris = connect.([(1, 5, 3), (4, 6, 2)], Triangle)
+    quads = connect.([(1, 2, 6, 5), (4, 3, 5, 6)], Quadrangle)
+    mesh = SimpleMesh(points, [tris; quads])
+
+    # Integrand & Solution
+    function integrand(p::Meshes.Point; a = a)
+        x₁, x₂ = ustrip.((to(p)))
+        (√(a^2 - x₁^2) + √(a^2 - x₂^2)) * u"A"
+    end
+    solution = 2a * (π * a^2 / 4) * u"A*m^2"
+
+    # Package and run tests
+    testable = TestableGeometry(integrand, mesh, solution)
+    runtests(testable; rtol = 1e-6)
+end
+
 @testitem "Meshes.Sphere 2D" setup=[Combinations] begin
     # Geometry
     origin = Point(0, 0)