diff --git a/CHANGELOG.md b/CHANGELOG.md index 0bc53834..69cecb3d 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/src/support.md b/docs/src/support.md index 0626a5fc..d6b57b9a 100644 --- a/docs/src/support.md +++ b/docs/src/support.md @@ -26,7 +26,7 @@ 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}` | 🛑 | ✅ | ✅ | @@ -34,6 +34,7 @@ designed to check for accuracy. | `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 | \ No newline at end of file +| 🛑 | Not supported | diff --git a/ext/MeshIntegralsEnzymeExt.jl b/ext/MeshIntegralsEnzymeExt.jl index 1b256c08..f78ad6b3 100644 --- a/ext/MeshIntegralsEnzymeExt.jl +++ b/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 diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 82da7fc4..201870d8 100644 --- a/src/MeshIntegrals.jl +++ b/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") diff --git a/src/integral.jl b/src/integral.jl index 1365aa93..94faebc5 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -3,7 +3,7 @@ ################################################################################ """ - 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 @@ -11,27 +11,44 @@ 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 ################################################################################ diff --git a/src/integral_aliases.jl b/src/integral_aliases.jl index 88e01bf3..811699ac 100644 --- a/src/integral_aliases.jl +++ b/src/integral_aliases.jl @@ -9,6 +9,9 @@ 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) @@ -16,18 +19,16 @@ Rule types available: """ 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,6 +42,9 @@ 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) @@ -48,18 +52,16 @@ Algorithm types available: """ 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,6 +75,9 @@ 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) @@ -80,16 +85,14 @@ Algorithm types available: """ 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 diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 5a8778d5..2de49fb3 100644 --- a/src/specializations/BezierCurve.jl +++ b/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. diff --git a/src/specializations/PolyArea.jl b/src/specializations/PolyArea.jl new file mode 100644 index 00000000..f2f93eff --- /dev/null +++ b/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 diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl index 12d21348..eb7a56fd 100644 --- a/src/specializations/Ring.jl +++ b/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 diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl index e917708b..9e6c536a 100644 --- a/src/specializations/Rope.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index 2b8d943c..94af3d2f 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/combinations.jl b/test/combinations.jl index d22a77e8..361b959e 100644 --- a/test/combinations.jl +++ b/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 + using MeshIntegrals: MeshIntegrals, GeometryOrDomain + 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 <: 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 <: 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) @@ -572,6 +592,25 @@ end runtests(testable) end +@testitem "Meshes.PolyArea" setup=[Combinations] begin + # Geometry + a, b, c = 0.4, 0.6, 1 + outer = [(0, 0), (c, 0), (c, c), (0, c)] + hole = [(a, a), (a, b), (b, b), (b, a)] + area = PolyArea([outer, hole]) + + # Integrand & Solution + function integrand(p::Meshes.Point) + x, y = ustrip.(u"m", to(p)) + 2x * u"A" + end + solution = (c^2 - (b - a) * (b^2 - a^2)) * u"A*m^2" + + # Package and run tests + testable = TestableGeometry(integrand, area, solution) + runtests(testable) +end + @testitem "Meshes.Pyramid" setup=[Combinations] begin if pkgversion(Meshes) >= v"0.52.12" # Geometry @@ -628,6 +667,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 +736,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 +747,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) @@ -727,6 +806,27 @@ end runtests(testable) end +@testitem "Meshes.StructuredGrid" setup=[Combinations] begin + # Geometry + a = π + N = 3 + side = range(0, a, length = N) + X = repeat(side, 1, N) + Y = repeat(side', N, 1) + grid = StructuredGrid(X, Y) + + # 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.Tetrahedron" setup=[Combinations] begin # Geometry a = Point(0, 0, 0)