v0.14-DEV: Refactor changing FP to a keyword argument (#91)

* Make uses of Line, Plane, and Ray explicit

* Tidying up

* Code cleanup

* Bugfix

* Bugfix

* Convert FP to kwarg with pass-through to workers

* Convert FP to kwarg

* Convert FP to kwarg

* Qualify use of Meshes.segments [skip ci]

* Convert FP to kwarg

* Convert FP to kwarg, wrapper no longer required for lineintegral

* Change comma to semicolon

* Convert FP to kwarg

* Update tests

* Fix erroneous leftover code fragment

* Apply suggestions from JuliaFormatter code review [skip ci]

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Tighten up line lengths

* Add missing kwargs

* Drop unbound type parameter

* apply format suggestions

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* default integration rules for `lineintegral`, `surfaceintegral`, and `volumeintegral`

* drop unbound type parameter

* fix default integration rule

* fix Floats to Vectors

* Qualify usages of Meshes area and volume

* Qualify usage of Meshes degree and clarify plane

* Remove unused explicit import - probably mistaken for I from IntegrationRule

* Convert Meshes and Unitful to imports

* Kill unitdirection function

* Split differentiation functions into separate file

* Improve docstring for jacobian

* Rename GaussKronrod-specific workers

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Joshua Lampert <[email protected]>
Co-authored-by: Joshua Lampert <[email protected]>

Author: 3 people

src/MeshIntegrals.jl

@@ -1,14 +1,16 @@
 module MeshIntegrals
 using CoordRefSystems: CoordRefSystems, CRS
-using LinearAlgebra: LinearAlgebra, I, norm, ×, ⋅
-using Meshes: Meshes, Line, Plane, Ray, area, degree, plane, segments, volume
-using Unitful: Unitful
+using LinearAlgebra: LinearAlgebra, norm, ×, ⋅
 
 import FastGaussQuadrature
 import HCubature
+import Meshes
 import QuadGK
+import Unitful
 
 include("utils.jl")
+
+include("differentiation.jl")
 export jacobian
 
 include("integration_rules.jl")

src/differentiation.jl

@@ -0,0 +1,125 @@
+################################################################################
+#                     Jacobian and Differential Elements
+################################################################################
+
+"""
+    jacobian(geometry, ts; ε=1e-6)
+
+Calculate the Jacobian of a geometry at some parametric point `ts` using a simple
+finite-difference approximation with step size `ε`.
+
+# Arguments
+- `geometry`: some `Meshes.Geometry` of N parametric dimensions
+- `ts`: a parametric point specified as a vector or tuple of length N
+- `ε`: step size to use for the finite-difference approximation
+"""
+function jacobian(
+        geometry,
+        ts;
+        ε = 1e-6
+)
+    T = eltype(ts)
+
+    # Get the partial derivative along the n'th axis via finite difference approximation
+    #   where ts is the current parametric position (εv is a reusable buffer)
+    function ∂ₙr!(εv, ts, n)
+        if ts[n] < T(0.01)
+            return ∂ₙr_right!(εv, ts, n)
+        elseif T(0.99) < ts[n]
+            return ∂ₙr_left!(εv, ts, n)
+        else
+            return ∂ₙr_central!(εv, ts, n)
+        end
+    end
+
+    # Central finite difference
+    function ∂ₙr_central!(εv, ts, n)
+        εv[n] = ε
+        a = ts .- εv
+        b = ts .+ εv
+        return (geometry(b...) - geometry(a...)) / 2ε
+    end
+
+    # Left finite difference
+    function ∂ₙr_left!(εv, ts, n)
+        εv[n] = ε
+        a = ts .- εv
+        b = ts
+        return (geometry(b...) - geometry(a...)) / ε
+    end
+
+    # Right finite difference
+    function ∂ₙr_right!(εv, ts, n)
+        εv[n] = ε
+        a = ts
+        b = ts .+ εv
+        return (geometry(b...) - geometry(a...)) / ε
+    end
+
+    # Allocate a re-usable ε vector
+    εv = zeros(T, length(ts))
+
+    ∂ₙr(n) = ∂ₙr!(εv, ts, n)
+    return map(∂ₙr, 1:length(ts))
+end
+
+"""
+    differential(geometry, ts)
+
+Calculate the differential element (length, area, volume, etc) of the parametric
+function for `geometry` at arguments `ts`.
+"""
+function differential(
+        geometry,
+        ts
+)
+    J = jacobian(geometry, ts)
+
+    # TODO generalize this with geometric algebra, e.g.: norm(foldl(∧, J))
+    if length(J) == 1
+        return norm(J[1])
+    elseif length(J) == 2
+        return norm(J[1] × J[2])
+    elseif length(J) == 3
+        return abs((J[1] × J[2]) ⋅ J[3])
+    else
+        error("Not implemented yet. Please report this as an Issue on GitHub.")
+    end
+end
+
+################################################################################
+#                             Partial Derivatives
+################################################################################
+
+"""
+    derivative(b::BezierCurve, t)
+
+Determine the vector derivative of a Bezier curve `b` for the point on the
+curve parameterized by value `t`.
+"""
+function derivative(
+        bz::Meshes.BezierCurve,
+        t
+)
+    # 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])
+    return N .* sum(sigma, 0:(N - 1))
+end

src/integral.jl

@@ -3,9 +3,7 @@
 ################################################################################
 
 """
-    integral(f, geometry)
-    integral(f, geometry, rule)
-    integral(f, geometry, rule, FP)
+    integral(f, geometry[, rule]; FP=Float64)
 
 Numerically integrate a given function `f(::Point)` over the domain defined by
 a `geometry` using a particular numerical integration `rule` with floating point
@@ -24,23 +22,13 @@ contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision
 function integral end
 
 # If only f and geometry are specified, select default rule
-function integral(
-        f::F,
-        geometry::G
-) where {F <: Function, G <: Meshes.Geometry}
-    N = Meshes.paramdim(geometry)
-    rule = (N == 1) ? GaussKronrod() : HAdaptiveCubature()
-    _integral(f, geometry, rule)
-end
-
-# with rule and T specified
 function integral(
         f::F,
         geometry::G,
-        rule::I,
-        FP::Type{T} = Float64
-) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule, T <: AbstractFloat}
-    _integral(f, geometry, rule, FP)
+        rule::I = Meshes.paramdim(geometry) == 1 ? GaussKronrod() : HAdaptiveCubature();
+        kwargs...
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule}
+    _integral(f, geometry, rule; kwargs...)
 end
 
 ################################################################################
@@ -51,25 +39,25 @@ end
 function _integral(
         f,
         geometry,
-        rule::GaussKronrod,
-        FP::Type{T} = Float64
-) where {T <: AbstractFloat}
-    # Run the appropriate integral type
+        rule::GaussKronrod;
+        kwargs...
+)
+    # Pass through to dim-specific workers in next section of this file
     N = Meshes.paramdim(geometry)
     if N == 1
-        return _integral_1d(f, geometry, rule, FP)
+        return _integral_gk_1d(f, geometry, rule; kwargs...)
     elseif N == 2
-        return _integral_2d(f, geometry, rule, FP)
+        return _integral_gk_2d(f, geometry, rule; kwargs...)
     elseif N == 3
-        return _integral_3d(f, geometry, rule, FP)
+        return _integral_gk_3d(f, geometry, rule; kwargs...)
     end
 end
 
 # GaussLegendre
 function _integral(
         f,
         geometry,
-        rule::GaussLegendre,
+        rule::GaussLegendre;
         FP::Type{T} = Float64
 ) where {T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
@@ -94,7 +82,7 @@ end
 function _integral(
         f,
         geometry,
-        rule::HAdaptiveCubature,
+        rule::HAdaptiveCubature;
         FP::Type{T} = Float64
 ) where {T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
@@ -103,7 +91,7 @@ function _integral(
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = fill(FP(0.5), 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(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
@@ -118,20 +106,20 @@ end
 #                    Specialized GaussKronrod Methods
 ################################################################################
 
-function _integral_1d(
+function _integral_gk_1d(
         f,
         geometry,
-        rule::GaussKronrod,
+        rule::GaussKronrod;
         FP::Type{T} = Float64
 ) where {T <: AbstractFloat}
     integrand(t) = f(geometry(t)) * differential(geometry, (t))
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
-function _integral_2d(
+function _integral_gk_2d(
         f,
         geometry2d,
-        rule::GaussKronrod,
+        rule::GaussKronrod;
         FP::Type{T} = Float64
 ) where {T <: AbstractFloat}
     integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v))
@@ -140,10 +128,10 @@ function _integral_2d(
 end
 
 # Integrating volumes with GaussKronrod not supported by default
-function _integral_3d(
+function _integral_gk_3d(
         f,
         geometry,
-        rule::GaussKronrod,
+        rule::GaussKronrod;
         FP::Type{T} = Float64
 ) where {T <: AbstractFloat}
     error("Integrating this volume type with GaussKronrod not supported.")