diff --git a/benchmark/Project.toml b/benchmark/Project.toml index 7ab23bea..96f687ac 100644 --- a/benchmark/Project.toml +++ b/benchmark/Project.toml @@ -2,8 +2,11 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] BenchmarkTools = "1.5" LinearAlgebra = "1" Meshes = "0.50, 0.51, 0.52" +Unitful = "1.19" +julia = "1.9" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index d830f578..f8f70ee9 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -2,6 +2,7 @@ using BenchmarkTools using LinearAlgebra using Meshes using MeshIntegrals +using Unitful const SUITE = BenchmarkGroup() @@ -19,10 +20,8 @@ rules = ( (name = "HAdaptiveCubature", rule = HAdaptiveCubature(), N = 500) ) geometries = ( - (name = "Meshes.BezierCurve", - item = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)])), - (name = "Meshes.Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))), - (name = "Meshes.Sphere", item = Sphere(Point(0, 0, 0), 1.0)) + (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))), + (name = "Sphere", item = Sphere(Point(0, 0, 0), 1.0)) ) SUITE["Integrals"] = let s = BenchmarkGroup() @@ -33,6 +32,43 @@ SUITE["Integrals"] = let s = BenchmarkGroup() s end +############################################################################################ +# Specializations +############################################################################################ + +spec = ( + f = p -> norm(to(p)), + f_exp = p::Point -> exp(-norm(to(p))^2 / u"m^2"), + g = ( + bezier = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)]), + line = Line(Point(0, 0, 0), Point(1, 1, 1)), + plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)), + ray = Ray(Point(0, 0, 0), Vec(0, 0, 1)), + triangle = Triangle(Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)), + tetrahedron = let + a = Point(0, 3, 0) + b = Point(-7, 0, 0) + c = Point(8, 0, 0) + ẑ = Vec(0, 0, 1) + Tetrahedron(a, b, c, a + ẑ) + end + ), + rule_gl = GaussLegendre(100), + rule_gk = GaussKronrod(), + rule_hc = HAdaptiveCubature() +) + +SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup() + s["BezierCurve"] = @benchmarkable integral($spec.f, $spec.g.bezier, $spec.rule_gl) + s["Line"] = @benchmarkable integral($spec.f_exp, $spec.g.line, $spec.rule_gl) + s["Plane"] = @benchmarkable integral($spec.f_exp, $spec.g.plane, $spec.rule_gl) + s["Ray"] = @benchmarkable integral($spec.f_exp, $spec.g.ray, $spec.rule_gl) + s["Triangle"] = @benchmarkable integral($spec.f, $spec.g.triangle, $spec.rule_gl) + # s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule) + # TODO re-enable Tetrahedron-GaussLegendre test when supported by main branch + s +end + ############################################################################################ # Differentials ############################################################################################ diff --git a/docs/src/supportmatrix.md b/docs/src/supportmatrix.md index 4265c06a..69234243 100644 --- a/docs/src/supportmatrix.md +++ b/docs/src/supportmatrix.md @@ -55,7 +55,7 @@ using an equivalent H-Adaptive Cubature rule, so support for this usage is not r | `SimpleMesh` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | | `Sphere` in `𝔼{2}` | ✅ | ✅ | ✅ | | `Sphere` in `𝔼{3}` | ✅ | ⚠️ | ✅ | -| `Tetrahedron` in `𝔼{3}` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | ✅ | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | +| `Tetrahedron` | ✅ | ⚠️ | ✅ | | `Triangle` | ✅ | ✅ | ✅ | | `Torus` | ✅ | ⚠️ | ✅ | | `Wedge` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 5fbac40f..7b270c7b 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -24,6 +24,7 @@ include("integral_aliases.jl") export lineintegral, surfaceintegral, volumeintegral # Integration methods specialized for particular geometries +include("specializations/_ParametricGeometry.jl") include("specializations/BezierCurve.jl") include("specializations/ConeSurface.jl") include("specializations/CylinderSurface.jl") diff --git a/src/integral.jl b/src/integral.jl index 4bb8ffc3..8777b86b 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -104,7 +104,7 @@ function _integral( # HCubature doesn't support functions that output Unitful Quantity types # Establish the units that are output by f - testpoint_parametriccoord = zeros(FP, N) + testpoint_parametriccoord = _zeros(FP, N) integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) # Create a wrapper that returns only the value component in those units uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts)) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index f8082e1b..01b0f24d 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -34,93 +34,22 @@ steep performance cost, especially for curves with a large number of control poi function integral( f, curve::Meshes.BezierCurve, - rule::GaussLegendre; - diff_method::DM = _default_method(curve), - FP::Type{T} = Float64, - alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, rule.n) - - # Change of variables: x [-1,1] ↦ t [0,1] - t(x) = (1 // 2) * x + (1 // 2) - point(x) = curve(t(x), alg) - integrand(x) = f(point(x)) * differential(curve, (t(x),), diff_method) - - # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length] - return (1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs)) -end - -function integral( - f, - curve::Meshes.BezierCurve, - rule::GaussKronrod; - diff_method::DM = _default_method(curve), - FP::Type{T} = Float64, - alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - point(t) = curve(t, alg) - integrand(t) = f(point(t)) * differential(curve, (t,), diff_method) - return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] -end - -function integral( - f, - curve::Meshes.BezierCurve, - rule::HAdaptiveCubature; - diff_method::DM = _default_method(curve), - FP::Type{T} = Float64, - alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - point(t) = curve(t, alg) - integrand(ts) = f(point(only(ts))) * differential(curve, ts, diff_method) - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = zeros(FP, 1) - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, _zeros(FP, 1), _ones(FP, 1); rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits + rule::IntegrationRule; + alg::Meshes.BezierEvalMethod = Meshes.Horner(), + kwargs... +) + # Generate a _ParametricGeometry whose parametric function auto-applies the alg kwarg + paramfunction(t) = _parametric(curve, t, alg) + param_curve = _ParametricGeometry(paramfunction, Meshes.paramdim(curve)) + + # Integrate the _ParametricGeometry using the standard methods + return _integral(f, param_curve, rule; kwargs...) end ################################################################################ -# jacobian +# Parametric ################################################################################ -function jacobian( - bz::Meshes.BezierCurve, - ts::V, - diff_method::Analytical -) where {V <: Union{AbstractVector, Tuple}} - t = only(ts) - # Parameter t restricted to domain [0,1] by definition - if t < 0 || t > 1 - throw(DomainError(t, "b(t) is not defined for t outside [0, 1].")) - end - - # Aliases - P = bz.controls - N = Meshes.degree(bz) - - # Ensure that this implementation is tractible: limited by ability to calculate - # binomial(N, N/2) without overflow. It's possible to extend this range by - # converting N to a BigInt, but this results in always returning BigFloat's. - N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.") - - # Generator for Bernstein polynomial functions - B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i) - - # Derivative = N Σ_{i=0}^{N-1} sigma(i) - # P indices adjusted for Julia 1-based array indexing - sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1]) - derivative = N .* sum(sigma, 0:(N - 1)) - - return (derivative,) +function _parametric(curve::Meshes.BezierCurve, t, alg::Meshes.BezierEvalMethod) + return curve(t, alg) end - -_has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index c877397a..2f6fdece 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -12,44 +12,32 @@ function integral( f, tetrahedron::Meshes.Tetrahedron, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _error_unsupported_combination("Tetrahedron", "GaussLegendre") -end - -function integral( - f, - tetrahedron::Meshes.Tetrahedron, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Tetrahedron, diff_method) - - o = zero(FP) - ∫uvw(u, v, w) = f(tetrahedron(u, v, w)) - ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), o, FP(1 - v - w); rule.kwargs...)[1] - ∫w(w) = QuadGK.quadgk(v -> ∫vw(v, w), o, FP(1 - w); rule.kwargs...)[1] - ∫ = QuadGK.quadgk(∫w, o, one(FP); rule.kwargs...)[1] - - # Apply barycentric domain correction (volume: 1/6 → actual) - return 6 * Meshes.volume(tetrahedron) * ∫ -end + rule::IntegrationRule; + kwargs... +) + # Generate a _ParametricGeometry whose parametric function domain spans [0,1]^3 + paramfunction(t1, t2, t3) = _parametric(tetrahedron, t1, t2, t3) + tetra = _ParametricGeometry(paramfunction, Meshes.paramdim(tetrahedron)) -function integral( - f, - tetrahedron::Meshes.Tetrahedron, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature") + # Integrate the _ParametricGeometry using the standard methods + return _integral(f, tetra, rule; kwargs...) end ################################################################################ -# jacobian +# Parametric ################################################################################ -_has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true +function _parametric(tetrahedron::Meshes.Tetrahedron, t1, t2, t3) + if (t3 < 0) || (t3 > 1) + msg = "tetrahedron(t1, t2, t3) is not defined for t3 outside [0, 1]." + throw(DomainError(t3, msg)) + end + + # Take a triangular cross-section at t3 + a = tetrahedron(t3, 0, 0) + b = tetrahedron(0, t3, 0) + c = tetrahedron(0, 0, t3) + cross_section = Meshes.Triangle(a, b, c) + + return _parametric(cross_section, t1, t2) +end diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index 28af5290..2c70f400 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -87,10 +87,18 @@ function integral( v = R * (1 - b / (a + b)) return f(triangle(u, v)) * R / (a + b)^2 end - ∫ = HCubature.hcubature(integrand, _zeros(FP, 2), (FP(1), FP(π / 2)), rule.kwargs...)[1] + + # HCubature doesn't support functions that output Unitful Quantity types + # Establish the units that are output by f + testpoint_parametriccoord = _zeros(T, 2) + integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) + # Create a wrapper that returns only the value component in those units + uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts)) + # Integrate only the unitless values + ∫ = HCubature.hcubature(uintegrand, _zeros(FP, 2), (T(1), T(π / 2)), rule.kwargs...)[1] # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * Meshes.area(triangle) .* ∫ + return 2 * Meshes.area(triangle) * integrandunits .* ∫ end ################################################################################ @@ -98,3 +106,21 @@ end ################################################################################ _has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true + +################################################################################ +# Parametric +################################################################################ + +function _parametric(triangle::Meshes.Triangle, t1, t2) + if (t1 < 0 || t1 > 1) || (t2 < 0 || t2 > 1) + msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]^2." + throw(DomainError((t1, t2), msg)) + end + + # Form a line segment between sides + a = triangle(0, t2) + b = triangle(t2, 0) + segment = Meshes.Segment(a, b) + + return segment(t1) +end diff --git a/src/specializations/_ParametricGeometry.jl b/src/specializations/_ParametricGeometry.jl new file mode 100644 index 00000000..16904d27 --- /dev/null +++ b/src/specializations/_ParametricGeometry.jl @@ -0,0 +1,63 @@ +""" + _ParametricGeometry <: Meshes.Primitive <: Meshes.Geometry + +_ParametricGeometry is used internally in MeshIntegrals.jl to behave like a generic wrapper +for geometries with custom parametric functions. This type is used for transforming other +geometries to enable integration over the standard rectangular `[0,1]^n` domain. + +Meshes.jl adopted a `ParametrizedCurve` type that performs a similar role as of `v0.51.20`, +but only supports geometries with one parametric dimension. Support is additionally planned +for more types that span surfaces and volumes, at which time this custom type will probably +no longer be required. + +# Fields +- `fun` - anything callable representing a parametric function: (ts...) -> Meshes.Point + +# Type Structure +- `M <: Meshes.Manifold` - same usage as in `Meshes.Geometry{M, C}` +- `C <: CoordRefSystems.CRS` - same usage as in `Meshes.Geometry{M, C}` +- `F` - type of the callable parametric function +- `Dim` - number of parametric dimensions +""" +struct _ParametricGeometry{M <: Meshes.Manifold, C <: CRS, F, Dim} <: + Meshes.Primitive{M, C} + fun::F + + function _ParametricGeometry{M, C}( + fun::F, + Dim::Int64 + ) where {M <: Meshes.Manifold, C <: CRS, F} + return new{M, C, F, Dim}(fun) + end +end + +""" + _ParametricGeometry(fun, dims) + +Construct a `_ParametricGeometry` using a provided parametric function `fun` for a geometry +with `dims` parametric dimensions. + +# Arguments +- `fun` - anything callable representing a parametric function mapping `(ts...) -> Meshes.Point` +- `dims::Int64` - number of parametric dimensions, i.e. `length(ts)` +""" +function _ParametricGeometry( + fun, + Dim::Int64 +) + p = fun(_zeros(Dim)...) + return _ParametricGeometry{Meshes.manifold(p), Meshes.crs(p)}(fun, Dim) +end + +# Allow a _ParametricGeometry to be called like a Geometry +(g::_ParametricGeometry)(ts...) = g.fun(ts...) + +Meshes.paramdim(::_ParametricGeometry{M, C, F, Dim}) where {M, C, F, Dim} = Dim + +""" + _parametric(geometry::G, ts...) where {G <: Meshes.Geometry} + +Used in MeshIntegrals.jl for defining parametric functions that transform non-standard +geometries into a form that can be integrated over the standard rectangular [0,1]^n limits. +""" +function _parametric end diff --git a/test/combinations.jl b/test/combinations.jl index 4dbd16dd..ed32a0a7 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -832,32 +832,32 @@ end @test_throws "not supported" volumeintegral(f, sphere) end -@testitem "Meshes.Tetrahedron" tags=[:extended] setup=[Setup] begin - pt_n = Point(0, 1, 0) - pt_w = Point(-1, 0, 0) - pt_e = Point(1, 0, 0) +@testitem "Meshes.Tetrahedron" setup=[Setup] begin + pt_n = Point(0, 3, 0) + pt_w = Point(-7, 0, 0) + pt_e = Point(8, 0, 0) ẑ = Vec(0, 0, 1) tetrahedron = Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ) - f(p) = 1.0 + f(p) = 1.0u"A" fv(p) = fill(f(p), 3) # Scalar integrand - sol = Meshes.measure(tetrahedron) - @test_throws "not supported" integral(f, tetrahedron, GaussLegendre(100)) - @test integral(f, tetrahedron, GaussKronrod()) ≈ sol - @test_throws "not supported" integral(f, tetrahedron, HAdaptiveCubature()) + sol = Meshes.measure(tetrahedron) * u"A" + @test integral(f, tetrahedron, GaussLegendre(100)) ≈ sol + @test_throws "not supported" integral(f, tetrahedron, GaussKronrod()) + @test integral(f, tetrahedron, HAdaptiveCubature()) ≈ sol # Vector integrand vsol = fill(sol, 3) - @test_throws "not supported" integral(fv, tetrahedron, GaussLegendre(100))≈vsol - @test integral(fv, tetrahedron, GaussKronrod()) ≈ vsol - @test_throws "not supported" integral(fv, tetrahedron, HAdaptiveCubature())≈vsol + @test integral(fv, tetrahedron, GaussLegendre(100)) ≈ vsol + @test_throws "not supported" integral(fv, tetrahedron, GaussKronrod()) + @test integral(fv, tetrahedron, HAdaptiveCubature()) ≈ vsol # Integral aliases @test_throws "not supported" lineintegral(f, tetrahedron) @test_throws "not supported" surfaceintegral(f, tetrahedron) - @test volumeintegral(f, tetrahedron, GaussKronrod()) ≈ sol + @test volumeintegral(f, tetrahedron, HAdaptiveCubature()) ≈ sol end @testitem "Meshes.Torus" setup=[Setup] begin @@ -892,11 +892,11 @@ end pt_e = Point(1, 0, 0) triangle = Triangle(pt_e, pt_n, pt_w) - f(p) = 1.0 + f(p) = 1.0u"A" fv(p) = fill(f(p), 3) # Scalar integrand - sol = Meshes.measure(triangle) + sol = Meshes.measure(triangle) * u"A" @test integral(f, triangle, GaussLegendre(100)) ≈ sol @test integral(f, triangle, GaussKronrod()) ≈ sol @test integral(f, triangle, HAdaptiveCubature()) ≈ sol diff --git a/test/utils.jl b/test/utils.jl index a42ae24f..69a8fde6 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -23,30 +23,53 @@ end using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical # _has_analytical of instances - bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) - @test _has_analytical(bezier) == true + triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0)) + @test _has_analytical(triangle) == true sphere = Sphere(Point(0, 0, 0), 1.0) @test _has_analytical(sphere) == false + # _default_method + @test _default_method(Meshes.Triangle) isa Analytical + @test _default_method(triangle) isa Analytical + @test _default_method(Meshes.Sphere) isa FiniteDifference + @test _default_method(sphere) isa FiniteDifference + # _has_analytical of types - @test _has_analytical(Meshes.BezierCurve) == true + @test _has_analytical(Meshes.BezierCurve) == false @test _has_analytical(Meshes.Line) == true @test _has_analytical(Meshes.Plane) == true @test _has_analytical(Meshes.Ray) == true @test _has_analytical(Meshes.Sphere) == false - @test _has_analytical(Meshes.Tetrahedron) == true + @test _has_analytical(Meshes.Tetrahedron) == false @test _has_analytical(Meshes.Triangle) == true # _guarantee_analytical @test _guarantee_analytical(Meshes.Line, Analytical()) === nothing @test_throws "Analytical" _guarantee_analytical(Meshes.Line, FiniteDifference()) - # _default_method - @test _default_method(Meshes.BezierCurve) isa Analytical - @test _default_method(bezier) isa Analytical - @test _default_method(Meshes.Sphere) isa FiniteDifference - @test _default_method(sphere) isa FiniteDifference - # FiniteDifference @test FiniteDifference().ε ≈ 1e-6 end + +@testitem "_ParametricGeometry" setup=[Setup] begin + using MeshIntegrals: _parametric + + pt_n = Point(0, 3, 0) + pt_w = Point(-7, 0, 0) + pt_e = Point(8, 0, 0) + ẑ = Vec(0, 0, 1) + triangle = Triangle(pt_n, pt_w, pt_e) + tetrahedron = Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ) + + # Ensure error is thrown for an out-of-bounds coordinate + for ts in [(-0.5, 0.5), (0.5, -0.5), (1.5, 0.5), (0.5, 1.5)] + @test_throws "not defined" _parametric(triangle, ts...) + # Tetrahedron forwards t1 and t2 to _parametric(::Triangle, ts...) + @test_throws "not defined" _parametric(tetrahedron, ts..., 0) + end + + # Ensue error is thrown for an out-of-bounds third coordinate + for t3 in [-0.5, 1.5] + @test_throws "not defined" _parametric(tetrahedron, 0, 0, t3) + end +end