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/specializations.md b/docs/src/specializations.md index dc873601..91a82e2a 100644 --- a/docs/src/specializations.md +++ b/docs/src/specializations.md @@ -14,6 +14,9 @@ There are several notable exceptions to how Meshes.jl defines [parametric functi ## Triangle +!!! note + This derivation is now obsolete. + For a specified `Meshes.Triangle` surface with area $A$, let $u$ and $v$ be Barycentric coordinates that span the surface. ```math \int_\triangle f(\bar{r}) \, \text{d}A diff --git a/docs/src/supportmatrix.md b/docs/src/supportmatrix.md index 82fcd5f3..eddae3c0 100644 --- a/docs/src/supportmatrix.md +++ b/docs/src/supportmatrix.md @@ -52,7 +52,7 @@ Cubature integration rules are recommended (and the default). | `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` in `𝔼{3}` | ✅ | 🛑 | ✅ | | `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 1a1c7878..f137859a 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -95,12 +95,12 @@ 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)) # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, zeros(FP, N), ones(FP, N); rule.kwargs...)[1] + value = HCubature.hcubature(uintegrand, _zeros(FP, N), _ones(FP, N); rule.kwargs...)[1] # Reapply units return value .* integrandunits diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index fc5905b3..9295dd6b 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -32,95 +32,21 @@ calculating Jacobians that are used to calculate differential elements steep performance cost, especially for curves with a large number of control points. """ function integral( - f::F, + f::Function, curve::Meshes.BezierCurve, - rule::GaussLegendre; - diff_method::DM = _default_method(curve), - FP::Type{T} = Float64, - alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {F <: Function, 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) = FP(1 // 2) * x + FP(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 FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs)) -end - -function integral( - f::F, - curve::Meshes.BezierCurve, - rule::GaussKronrod; - diff_method::DM = _default_method(curve), - FP::Type{T} = Float64, - alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {F <: Function, 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::F, - curve::Meshes.BezierCurve, - rule::HAdaptiveCubature; - diff_method::DM = _default_method(curve), - FP::Type{T} = Float64, - alg::Meshes.BezierEvalMethod = Meshes.Horner() -) where {F <: Function, 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... +) + paramfunction(t) = _parametric(curve, t; alg = alg) + param_curve = _ParametricGeometry(paramfunction, 1) + 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/Line.jl b/src/specializations/Line.jl index 4d4c4c89..97082f69 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -8,83 +8,22 @@ ################################################################################ function integral( - f::F, + f::Function, line::Meshes.Line, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Line, diff_method) - - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, rule.n) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Line - differential(line, x) = t′(x) * _units(line(0)) - integrand(x) = f(line(t(x))) * differential(line, x) - return sum(w .* integrand(x) for (w, x) in zip(ws, xs)) -end - -function integral( - f::F, - line::Meshes.Line, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Line, diff_method) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Integrate f along the Line - domainunits = _units(line(0)) - integrand(t) = f(line(t)) * domainunits - return QuadGK.quadgk(integrand, FP(-Inf), FP(Inf); rule.kwargs...)[1] -end - -function integral( - f::F, - line::Meshes.Line, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Line, diff_method) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Line - differential(line, x) = t′(x) * _units(line(0)) - integrand(x::AbstractVector) = f(line(t(x[1]))) * differential(line, x[1]) - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = FP[0.5] - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, (-one(FP),), (one(FP),); rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits + rule::IntegrationRule; + kwargs... +) + paramfunction(t) = _parametric(line, t) + param_line = _ParametricGeometry(paramfunction, 1) + return _integral(f, param_line, rule; kwargs...) end ################################################################################ -# jacobian +# Parametric ################################################################################ -_has_analytical(::Type{T}) where {T <: Meshes.Line} = true +function _parametric(line::Meshes.Line, t) + f1(t) = t / (1 - t^2) + f2(t) = 2t - 1 + return line((f1 ∘ f2)(t)) +end diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index 1b34405e..fd0b56d6 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -1,104 +1,30 @@ -################################################################################ -# Specialized Methods for Plane +############################################################################################ +# Specialized Methods for Plane # # Why Specialized? -# The Plane geometry is a special case, representing a planar surface with -# infinite extent along two basis vectors. This requires a pair of domain -# transformations mapping from the typical parametric region [0,1]² to an -# infinite one (-∞,∞)². -################################################################################ +# The Plane geometry is a special case, representing a planar surface with infinite extent +# along two basis vectors. This requires a pair of domain transformations mapping from the +# typical parametric region [0,1]² to an infinite one (-∞,∞)². +############################################################################################ function integral( - f::F, + f::Function, plane::Meshes.Plane, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Plane, diff_method) - - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² - xs, ws = _gausslegendre(FP, rule.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Normalize the Plane's orthogonal vectors - uu = Meshes.unormalize(plane.u) - uv = Meshes.unormalize(plane.v) - uplane = Meshes.Plane(plane.p, uu, uv) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f over the Plane - domainunits = _units(uplane(0, 0)) - function integrand(((wi, wj), (xi, xj))) - wi * wj * f(uplane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2 - end - return sum(integrand, zip(wws, xxs)) -end - -function integral( - f::F, - plane::Meshes.Plane, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Plane, diff_method) - - # Normalize the Plane's orthogonal vectors - uu = Meshes.unormalize(plane.u) - uv = Meshes.unormalize(plane.v) - uplane = Meshes.Plane(plane.p, uu, uv) - - # Integrate f over the Plane - domainunits = _units(uplane(0, 0))^2 - integrand(u, v) = f(uplane(u, v)) * domainunits - inner∫(v) = QuadGK.quadgk(u -> integrand(u, v), FP(-Inf), FP(Inf); rule.kwargs...)[1] - return QuadGK.quadgk(inner∫, FP(-Inf), FP(Inf); rule.kwargs...)[1] + rule::IntegrationRule; + kwargs... +) + paramfunction(t1, t2) = _parametric(plane, t1, t2) + param_plane = _ParametricGeometry(paramfunction, 2) + return _integral(f, param_plane, rule; kwargs...) end -function integral( - f::F, - plane::Meshes.Plane, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Plane, diff_method) - - # Normalize the Plane's orthogonal vectors - uu = Meshes.unormalize(plane.u) - uv = Meshes.unormalize(plane.v) - uplane = Meshes.Plane(plane.p, uu, uv) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 +############################################################################################ +# Parametric +############################################################################################ - # Integrate f over the Plane - domainunits = _units(uplane(0, 0)) - function integrand(x::AbstractVector) - f(uplane(t(x[1]), t(x[2]))) * t′(x[1]) * t′(x[2]) * domainunits^2 - end - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = zeros(FP, 2) - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, -ones(FP, 2), ones(FP, 2); rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits +function _parametric(plane::Meshes.Plane, t1, t2) + f1(t) = t / (1 - t^2) + f2(t) = 2t - 1 + f = f1 ∘ f2 + return plane(f(t1), f(t2)) end - -################################################################################ -# jacobian -################################################################################ - -_has_analytical(::Type{T}) where {T <: Meshes.Plane} = true diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index 9008b4a8..67d0234c 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -1,95 +1,28 @@ -################################################################################ -# Specialized Methods for Ray +############################################################################################ +# Specialized Methods for Ray # # Why Specialized? -# The Ray geometry is a special case, representing a ray, originating at a point -# and extending an infinite length in a particular direction. This requires -# a domain transformation mapping from the typical parametric region [0,1] to -# an infinite one (-∞,∞). -################################################################################ +# The Ray geometry is a special case, representing a ray, originating at a point and +# extending an infinite length in a particular direction. This requires a domain transform +# that maps from the typical parametric region [0,1] to an infinite one (-∞,∞). +############################################################################################ function integral( - f::F, + f::Function, ray::Meshes.Ray, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Ray, diff_method) - - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, rule.n) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞) - t₁(x) = FP(1 / 2) * x + FP(1 / 2) - t₁′(x) = FP(1 / 2) - t₂(x) = x / (1 - x^2) - t₂′(x) = (1 + x^2) / (1 - x^2)^2 - t = t₂ ∘ t₁ - t′(x) = t₂′(t₁(x)) * t₁′(x) - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(x) = f(ray(t(x))) * t′(x) * domainunits - return sum(w .* integrand(x) for (w, x) in zip(ws, xs)) -end - -function integral( - f::F, - ray::Meshes.Ray, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Ray, diff_method) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(t) = f(ray(t)) * domainunits - return QuadGK.quadgk(integrand, zero(FP), FP(Inf); rule.kwargs...)[1] + rule::IntegrationRule; + kwargs... +) + paramfunction(t) = _parametric(ray, t) + param_ray = _ParametricGeometry(paramfunction, 1) + return _integral(f, param_ray, rule; kwargs...) end -function integral( - f::F, - ray::Meshes.Ray, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Ray, diff_method) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) - - # Domain transformation: x ∈ [0,1] ↦ t ∈ [0,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 +############################################################################################ +# Parametric +############################################################################################ - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(x::AbstractVector) = f(ray(t(x[1]))) * t′(x[1]) * domainunits - - # 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(uv) = Unitful.ustrip.(integrandunits, integrand(uv)) - # Integrate only the unitless values - value = HCubature.hcubature(uintegrand, zeros(FP, 1), ones(FP, 1); rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits +function _parametric(ray::Meshes.Ray, t) + f(t) = t / (1 - t^2) + return ray(f(t)) end - -################################################################################ -# jacobian -################################################################################ - -_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index 9fb2af9a..d777f499 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -1,55 +1,39 @@ -################################################################################ +############################################################################################ # Specialized Methods for Tetrahedron # # Why Specialized? -# The Tetrahedron geometry is a volumetric simplex whose parametric function -# in Meshes.jl uses barycentric coordinates on a domain {u,v,w} with coordinates -# that are non-negative and bound by the surface $u + v + w ≤ 1$. This requires -# a multi-step domain transformation whose derivation is detailed in the package -# documentation. -################################################################################ +# The Tetrahedron geometry is a volumetric simplex whose parametric function in Meshes.jl +# uses barycentric coordinates on a domain {u,v,w} with coordinates that are non-negative +# and bounded by the surface $u + v + w ≤ 1$. A transformation is used to map this volume +# to a rectangular domain [0,1]^3. +############################################################################################ function integral( f::F, tetrahedron::Meshes.Tetrahedron, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _error_unsupported_combination("Tetrahedron", "GaussLegendre") + rule::IntegrationRule; + kwargs... +) where {F <: Function} + paramfunction(t1, t2, t3) = _parametric(tetrahedron, t1, t2, t3) + tetra = _ParametricGeometry(paramfunction, 3) + return _integral(f, tetra, rule; kwargs...) end -function integral( - f::F, - tetrahedron::Meshes.Tetrahedron, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Tetrahedron, diff_method) +################################################################################ +# Parametric +################################################################################ - nil = zero(FP) - ∫uvw(u, v, w) = f(tetrahedron(u, v, w)) - ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), nil, FP(1 - v - w); rule.kwargs...)[1] - ∫w(w) = QuadGK.quadgk(v -> ∫vw(v, w), nil, FP(1 - w); rule.kwargs...)[1] - ∫ = QuadGK.quadgk(∫w, nil, one(FP); rule.kwargs...)[1] +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 - # Apply barycentric domain correction (volume: 1/6 → actual) - return 6 * Meshes.volume(tetrahedron) * ∫ -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) -function integral( - f::F, - tetrahedron::Meshes.Tetrahedron, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature") + return _parametric(cross_section, t1, t2) end - -################################################################################ -# jacobian -################################################################################ - -_has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index c63e715b..7e9407a6 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -1,100 +1,38 @@ -################################################################################ +############################################################################################ # Specialized Methods for Triangle # # Why Specialized? -# The Triangle geometry is a surface simplex whose parametric function in -# Meshes.jl uses barycentric coordinates on a domain {u,v} with coordinates -# that are non-negative and bound by the surface $u + v ≤ 1$. This requires a -# multi-step domain transformation whose derivation is detailed in the package -# documentation. -################################################################################ - -function integral( - f::F, - triangle::Meshes.Triangle, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Triangle, diff_method) - - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 - xs, ws = _gausslegendre(FP, rule.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Domain transformations: - # xᵢ [-1,1] ↦ R [0,1] - # xⱼ [-1,1] ↦ φ [0,π/2] - uR(xᵢ) = T(1 // 2) * (xᵢ + 1) - uφ(xⱼ) = T(π / 4) * (xⱼ + 1) - - # Integrate the Barycentric triangle by transforming it into polar coordinates - # with a modified radius - # R = r ( sinφ + cosφ ) - # s.t. integration bounds become rectangular - # R ∈ [0, 1] and φ ∈ [0, π/2] - function integrand(((wᵢ, wⱼ), (xᵢ, xⱼ))) - R = uR(xᵢ) - φ = uφ(xⱼ) - a, b = sincos(φ) - u = R * (1 - a / (a + b)) - v = R * (1 - b / (a + b)) - return wᵢ * wⱼ * f(triangle(u, v)) * R / (a + b)^2 - end - - # Calculate 2D Gauss-Legendre integral over modified-polar-Barycentric coordinates - # Apply a linear domain-correction factor - return FP(π / 4) * Meshes.area(triangle) .* sum(integrand, zip(wws, xxs)) -end +# The Triangle geometry is a surface simplex whose parametric function in Meshes.jl uses +# barycentric coordinates on a domain {u,v} with coordinates that are non-negative and +# bound by the surface $u + v ≤ 1$. A transformation is used to map this volume to a +# rectangular domain [0,1]^3. +############################################################################################ function integral( f::F, triangle::Meshes.Triangle, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Triangle, diff_method) - - # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0) - # i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du - ∫u(u) = QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1] - ∫ = QuadGK.quadgk(∫u, zero(FP), one(FP); rule.kwargs...)[1] - - # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * Meshes.area(triangle) .* ∫ + rule::IntegrationRule; + kwargs... +) where {F <: Function} + paramfunction(t1, t2) = _parametric(triangle, t1, t2) + tri = _ParametricGeometry(paramfunction, 2) + return _integral(f, tri, rule; kwargs...) end -function integral( - f::F, - triangle::Meshes.Triangle, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {F <: Function, DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Triangle, diff_method) +################################################################################ +# Parametric +################################################################################ - # Integrate the Barycentric triangle by transforming it into polar coordinates - # with a modified radius - # R = r ( sinφ + cosφ ) - # s.t. integration bounds become rectangular - # R ∈ [0, 1] and φ ∈ [0, π/2] - function integrand(Rφ) - R, φ = Rφ - a, b = sincos(φ) - u = R * (1 - a / (a + b)) - v = R * (1 - b / (a + b)) - return f(triangle(u, v)) * R / (a + b)^2 +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 - ∫ = HCubature.hcubature(integrand, zeros(FP, 2), FP[1, π / 2], rule.kwargs...)[1] - - # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * Meshes.area(triangle) .* ∫ -end -################################################################################ -# jacobian -################################################################################ + # Form a line segment between sides + a = triangle(0, t2) + b = triangle(t2, 0) + segment = Meshes.Segment(a, b) -_has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true + return segment(t1) +end diff --git a/src/specializations/_ParametricGeometry.jl b/src/specializations/_ParametricGeometry.jl new file mode 100644 index 00000000..b5e8fadc --- /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::Function` - 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 integrand 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::Function` - parametric function mapping `(ts...) -> Meshes.Point` +- `dims::Int64` - number of parametric dimensions, i.e. `length(ts)` +""" +function _ParametricGeometry( + fun::F, + Dim::Int64 +) where {F <: Function} + 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/src/utils.jl b/src/utils.jl index e98ade74..c1cbc76f 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,20 +14,18 @@ function _error_unsupported_combination(geometry, rule) throw(ArgumentError(msg)) end +# Return an NTuple{N, T} of zeros; same interface as Base.zeros() but faster +_zeros(T::DataType, N::Int64) = ntuple(_ -> zero(T), N) +_zeros(N::Int) = _zeros(Float64, N) + +# Return an NTuple{N, T} of ones; same interface as Base.ones() but faster +_ones(T::DataType, N::Int64) = ntuple(_ -> one(T), N) +_ones(N::Int) = _ones(Float64, N) + ################################################################################ # DifferentiationMethod ################################################################################ -# Throw an ArgumentError if Analytical() jacobian not defined for this type -function _guarantee_analytical( - ::Type{G}, - ::DifferentiationMethod -) where {G <: Geometry} - throw(ArgumentError("$G geometries require kwarg diff_method = Analytical()")) -end - -_guarantee_analytical(::Type{G}, ::Analytical) where {G <: Geometry} = nothing - # Return whether a geometry type has jacobian methods defined _has_analytical(::Type{G}) where {G <: Geometry} = false _has_analytical(g::G) where {G <: Geometry} = _has_analytical(G) diff --git a/test/combinations.jl b/test/combinations.jl index aa5e87f5..38d3a9e4 100644 --- a/test/combinations.jl +++ b/test/combinations.jl @@ -829,32 +829,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())≈sol + @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())≈vsol + @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 @@ -884,16 +884,16 @@ end end @testitem "Meshes.Triangle" setup=[Setup] begin - pt_n = Point(0, 1, 0) - pt_w = Point(-1, 0, 0) - pt_e = Point(1, 0, 0) - triangle = Triangle(pt_e, pt_n, pt_w) + a(T = Float64) = Point(0, T(1), 0) + b(T = Float64) = Point(T(-1), 0, 0) + c(T = Float64) = Point(T(1), 0, 0) + triangle = Triangle(a(), b(), c()) - 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 @@ -908,4 +908,10 @@ end @test_throws "not supported" lineintegral(f, triangle) @test surfaceintegral(f, triangle) ≈ sol @test_throws "not supported" volumeintegral(f, triangle) + + # Type stability + f32(p) = 1.0f0u"A" + triangle32 = Triangle(a(Float32), b(Float32), c(Float32)) + int_FP32 = integral(f32, triangle32, GaussLegendre(100); FP = Float32) + @test typeof(int_FP32.val) == Float32 end diff --git a/test/utils.jl b/test/utils.jl index f5c26308..b290917a 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,5 +1,6 @@ @testitem "Utilities" setup=[Setup] begin using LinearAlgebra: norm + using MeshIntegrals: _units, _zeros, _ones # _KVector v = Meshes.Vec(3, 4) @@ -7,37 +8,58 @@ # _units p = Point(1.0u"cm", 2.0u"mm", 3.0u"m") - @test MeshIntegrals._units(p) == u"m" + @test _units(p) == u"m" + + # _zeros + @test _zeros(2) == (0.0, 0.0) + @test _zeros(Float32, 2) == (0.0f0, 0.0f0) + + # _ones + @test _ones(2) == (1.0, 1.0) + @test _ones(Float32, 2) == (1.0f0, 1.0f0) end @testitem "DifferentiationMethod" setup=[Setup] begin - using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical + using MeshIntegrals: _has_analytical, _default_method # _has_analytical of instances bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)]) - @test _has_analytical(bezier) == true + @test _has_analytical(bezier) == false sphere = Sphere(Point(0, 0, 0), 1.0) @test _has_analytical(sphere) == false # _has_analytical of types - @test _has_analytical(Meshes.BezierCurve) == true - @test _has_analytical(Meshes.Line) == true - @test _has_analytical(Meshes.Plane) == true - @test _has_analytical(Meshes.Ray) == true + @test _has_analytical(Meshes.BezierCurve) == false @test _has_analytical(Meshes.Sphere) == false - @test _has_analytical(Meshes.Tetrahedron) == true - @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.BezierCurve) isa FiniteDifference + @test _default_method(bezier) isa FiniteDifference @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: _ParametricGeometry, _parametric + + # paramdim(::_ParametricGeometry) + segment = Segment(Point(0, 0), Point(1, 1)) + f(t) = segment(t) + geometry = _ParametricGeometry(f, 1) + @test paramdim(geometry) == 1 + + # _parametric bounds checks + triangle = Triangle(Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)) + @test_throws DomainError _parametric(triangle, 1.1, 0.0) + 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 + @test_throws DomainError _parametric(tetrahedron, 1.1, 0.0, 0.0) +end