Use a new _ParametricGeometry for specializations

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

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) |

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")

src/specializations/BezierCurve.jl

@@ -32,60 +32,23 @@ 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]
+        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
 
-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]
+################################################################################
+#                              Parametric
+################################################################################
 
-    # Reapply units
-    return value .* integrandunits
+function _parametric(curve::Meshes.BezierCurve, t; alg::Meshes.BezierEvalMethod)
+    return curve(t, alg)
 end
 
 ################################################################################
@@ -94,9 +57,9 @@ end
 
 function jacobian(
         bz::Meshes.BezierCurve,
-        ts::V,
+        ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
         diff_method::Analytical
-) where {V <: Union{AbstractVector, Tuple}}
+) where {T <: AbstractFloat}
     t = only(ts)
     # Parameter t restricted to domain [0,1] by definition
     if t < 0 || t > 1

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

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