Use recently-implemented Meshes predicates

Project.toml

@@ -25,7 +25,7 @@ Enzyme = "0.13.47"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 LinearAlgebra = "1"
-Meshes = "0.52.12, 0.53, 0.54"
+Meshes = "0.53, 0.54"
 QuadGK = "2.1.1"
 Unitful = "1.19"
 julia = "1.9"

benchmark/Project.toml

@@ -3,12 +3,14 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 BenchmarkTools = "1.5"
-Enzyme = "0.13.19"
+Enzyme = "0.13.47"
 LinearAlgebra = "1"
-Meshes = "0.51.20, 0.52"
+Meshes = "0.53, 0.54"
+StaticArrays = "1"
 Unitful = "1.19"
 julia = "1.9"

benchmark/benchmarks.jl

@@ -2,6 +2,7 @@ using BenchmarkTools
 using LinearAlgebra
 using Meshes
 using MeshIntegrals
+using StaticArrays
 using Unitful
 import Enzyme
 
@@ -11,9 +12,12 @@ const SUITE = BenchmarkGroup()
 #                                      Integrals
 ############################################################################################
 
+# Like fill but returns an SVector to eliminate integrand function allocations
+staticfill(val, N) = SVector(ntuple(Returns(val), N)...)
+
 integrands = (
     (name = "Scalar", f = p -> norm(to(p))),
-    (name = "Vector", f = p -> fill(norm(to(p)), 3))
+    (name = "Vector", f = p -> staticfill(norm(to(p)), 3))
 )
 rules = (
     (name = "GaussLegendre", rule = GaussLegendre(100)),

docs/Project.toml

@@ -8,6 +8,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 BenchmarkTools = "1"
 Documenter = "1"
-Meshes = "0.51.20, 0.52"
+Meshes = "0.53, 0.54"
 Unitful = "1.19"
 julia = "1.9"

src/integral.jl

@@ -2,6 +2,9 @@
 #                         Master Integral Function
 ################################################################################
 
+# Default integration rule when rule arg is unspecified
+_default_rule(g) = ifelse(Meshes.iscurve(g), GaussKronrod(), HAdaptiveCubature())
+
 """
     integral(f, geometry[, rule]; kwargs...)
 
@@ -22,12 +25,6 @@ calculate Jacobians within the integration domain.
 """
 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,
@@ -64,7 +61,7 @@ function _integral(
     _check_diff_method_support(geometry, diff_method)
 
     # Only supported for 1D geometries
-    if Meshes.paramdim(geometry) != 1
+    if !Meshes.iscurve(geometry)
         msg = "GaussKronrod rules not supported in more than one parametric dimension."
         throw(ArgumentError(msg))
     end

src/integral_aliases.jl

@@ -1,5 +1,5 @@
 function _alias_error_msg(name, N)
-    "Performing a $name integral on a geometry with $N parametric dimensions not supported."
+    "$name integrals not supported on a geometry without exactly $N parametric dimensions."
 end
 
 ################################################################################
@@ -27,8 +27,8 @@ function lineintegral(
         rule::IntegrationRule = GaussKronrod();
         kwargs...
 )
-    if (N = Meshes.paramdim(geometry)) != 1
-        throw(ArgumentError(_alias_error_msg("line", N)))
+    if !Meshes.iscurve(geometry)
+        throw(ArgumentError(_alias_error_msg("Line", 1)))
     end
 
     return integral(f, geometry, rule; kwargs...)
@@ -59,8 +59,8 @@ function surfaceintegral(
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
 )
-    if (N = Meshes.paramdim(geometry)) != 2
-        throw(ArgumentError(_alias_error_msg("surface", N)))
+    if !Meshes.issurface(geometry)
+        throw(ArgumentError(_alias_error_msg("Surface", 2)))
     end
 
     return integral(f, geometry, rule; kwargs...)
@@ -91,8 +91,8 @@ function volumeintegral(
         rule::IntegrationRule = HAdaptiveCubature();
         kwargs...
 )
-    if (N = Meshes.paramdim(geometry)) != 3
-        throw(ArgumentError(_alias_error_msg("volume", N)))
+    if !Meshes.issolid(geometry)
+        throw(ArgumentError(_alias_error_msg("Volume", 3)))
     end
 
     return integral(f, geometry, rule; kwargs...)

test/Project.toml

@@ -6,19 +6,21 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-Aqua = "0.7, 0.8"
-CoordRefSystems = "0.15, 0.16, 0.17, 0.18"
+Aqua = "0.8"
+CoordRefSystems = "0.16, 0.17, 0.18"
 Enzyme = "0.13.47"
 ExplicitImports = "1.6.0"
 LinearAlgebra = "1"
-Meshes = "0.51.20, 0.52, 0.53, 0.54"
+Meshes = "0.53, 0.54"
 SpecialFunctions = "2"
+StaticArrays = "1"
 TestItemRunner = "1"
 TestItems = "1"
 Unitful = "1.19"

test/combinations.jl

@@ -18,9 +18,12 @@ This file includes tests for:
     import LinearAlgebra: norm
     using Meshes
     using MeshIntegrals: MeshIntegrals, GeometryOrDomain
+    import StaticArrays: SVector
     import Unitful: @u_str, Quantity, ustrip
     import Enzyme
 
+    staticfill(val, N) = SVector(ntuple(Returns(val), N)...)
+
     # Used for testing callable objects as integrand functions
     struct Callable{F <: Function}
         f::F
@@ -113,8 +116,8 @@ This file includes tests for:
                 @test integral(f, testable.geometry, rule)≈sol rtol=rtol
 
                 # Vector integrand
-                fv(p) = fill(testable.integrand(p), 3)
-                sol_v = fill(testable.solution, 3)
+                fv(p) = staticfill(testable.integrand(p), 3)
+                sol_v = staticfill(testable.solution, 3)
                 @test integral(fv, testable.geometry, rule)≈sol_v rtol=rtol
             else
                 f = testable.integrand