MeshIntegrals.jl/src/differentiation.jl at 8105b716d3f24d788e7b4008cef1c6049580b0a4 · JuliaGeometry/MeshIntegrals.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
################################################################################
#                          DifferentiationMethods
################################################################################

"""
    DifferentiationMethod

A category of types used to specify the desired method for calculating derivatives.
Derivatives are used to form Jacobian matrices when calculating the differential
element size throughout the integration region.

See also [`FiniteDifference`](@ref), [`Analytical`](@ref).
"""
abstract type DifferentiationMethod end

"""
    FiniteDifference(ε=1e-6)

Use to specify use of a finite-difference approximation method with a step size
of `ε` for calculating derivatives.
"""
struct FiniteDifference{T <: AbstractFloat} <: DifferentiationMethod
    ε::T
end

# Default constructors
FiniteDifference{T}() where {T <: AbstractFloat} = FiniteDifference{T}(T(1e-6))
FiniteDifference() = FiniteDifference{Float64}()

"""
    Analytical()

Use to specify use of analytically-derived solutions for calculating derivatives.
These solutions are currently defined only for a subset of geometry types.

# Supported Geometries:
- `BezierCurve`
- `Line`
- `Plane`
- `Ray`
- `Tetrahedron`
- `Triangle`
"""
struct Analytical <: DifferentiationMethod end

# Future Support:
#   struct AutoEnzyme <: DifferentiationMethod end
#   struct AutoZygote <: DifferentiationMethod end

################################################################################
#                                  Jacobian
################################################################################

"""
    jacobian(geometry, ts[, diff_method])

Calculate the Jacobian of a `geometry`'s parametric function at some point `ts`.
Optionally, direct the use of a particular differentiation method `diff_method`; by default
use analytic solutions where possible and finite difference approximations
otherwise.

# Arguments
- `geometry`: some `Meshes.Geometry` of N parametric dimensions
- `ts`: a parametric point specified as a vector or tuple of length N
- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use
"""
function jacobian(
        geometry::G,
        ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}
) where {G <: Geometry, T <: AbstractFloat}
    return jacobian(geometry, ts, _default_method(G))
end

function jacobian(
        geometry::Geometry,
        ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
        diff_method::FiniteDifference
) where {T <: AbstractFloat}
    Dim = Meshes.paramdim(geometry)
    if Dim != length(ts)
        throw(ArgumentError("ts must have same number of dimensions as geometry."))
    end

    # Get the partial derivative along the n'th axis via finite difference
    #   approximation, where ts is the current parametric position
    function ∂ₙr(ts, n, ε)
        # Build left/right parametric coordinates with non-allocating iterators
        left = Iterators.map(((i, t),) -> i == n ? t - ε : t, enumerate(ts))
        right = Iterators.map(((i, t),) -> i == n ? t + ε : t, enumerate(ts))
        # Select orientation of finite-diff
        if ts[n] < 0.01
            # Right
            return (geometry(right...) - geometry(ts...)) / ε
        elseif 0.99 < ts[n]
            # Left
            return (geometry(ts...) - geometry(left...)) / ε
        else
            # Central
            return (geometry(right...) - geometry(left...)) / 2ε
        end
    end

    return ntuple(n -> ∂ₙr(ts, n, T(diff_method.ε)), Dim)
end

################################################################################
#                          Differential Elements
################################################################################

"""
    differential(geometry, ts[, diff_method])

Calculate the differential element (length, area, volume, etc) of the parametric
function for `geometry` at arguments `ts`. Optionally, direct the use of a
particular differentiation method `diff_method`; by default use analytic solutions where
possible and finite difference approximations otherwise.

# Arguments
- `geometry`: some `Meshes.Geometry` of N parametric dimensions
- `ts`: a parametric point specified as a vector or tuple of length N
- `diff_method`: the desired [`DifferentiationMethod`](@ref) to use
"""
function differential(
        geometry::G,
        ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
        diff_method::DifferentiationMethod = _default_method(G)
) where {G <: Geometry, T <: AbstractFloat}
    J = Iterators.map(_KVector, jacobian(geometry, ts, diff_method))
    return LinearAlgebra.norm(foldl(∧, J))
end