@@ -25,36 +25,35 @@ function integral end
25
25
26
26
# If only f and geometry are specified, select default rule
27
27
function integral (
28
- f:: F ,
29
- geometry:: G
30
- ) where {F<: Function , G<: Meshes.Geometry }
28
+ f:: F ,
29
+ geometry:: G
30
+ ) where {F <: Function , G <: Meshes.Geometry }
31
31
N = Meshes. paramdim (geometry)
32
32
rule = (N == 1 ) ? GaussKronrod () : HAdaptiveCubature ()
33
33
_integral (f, geometry, rule)
34
34
end
35
35
36
36
# with rule and T specified
37
37
function integral (
38
- f:: F ,
39
- geometry:: G ,
40
- rule:: I ,
41
- FP:: Type{T} = Float64
42
- ) where {F<: Function , G<: Meshes.Geometry , I<: IntegrationRule , T<: AbstractFloat }
38
+ f:: F ,
39
+ geometry:: G ,
40
+ rule:: I ,
41
+ FP:: Type{T} = Float64
42
+ ) where {F <: Function , G <: Meshes.Geometry , I <: IntegrationRule , T <: AbstractFloat }
43
43
_integral (f, geometry, rule, FP)
44
44
end
45
45
46
-
47
46
# ###############################################################################
48
47
# Generalized (n-Dimensional) Worker Methods
49
48
# ###############################################################################
50
49
51
50
# GaussKronrod
52
51
function _integral (
53
- f,
54
- geometry,
55
- rule:: GaussKronrod ,
56
- FP:: Type{T} = Float64
57
- ) where {T<: AbstractFloat }
52
+ f,
53
+ geometry,
54
+ rule:: GaussKronrod ,
55
+ FP:: Type{T} = Float64
56
+ ) where {T <: AbstractFloat }
58
57
# Run the appropriate integral type
59
58
N = Meshes. paramdim (geometry)
60
59
if N == 1
68
67
69
68
# GaussLegendre
70
69
function _integral (
71
- f,
72
- geometry,
73
- rule:: GaussLegendre ,
74
- FP:: Type{T} = Float64
75
- ) where {T<: AbstractFloat }
70
+ f,
71
+ geometry,
72
+ rule:: GaussLegendre ,
73
+ FP:: Type{T} = Float64
74
+ ) where {T <: AbstractFloat }
76
75
N = Meshes. paramdim (geometry)
77
76
78
77
# Get Gauss-Legendre nodes and weights for a region [-1,1]^N
@@ -81,23 +80,23 @@ function _integral(
81
80
nodes = Iterators. product (ntuple (Returns (xs), N)... )
82
81
83
82
# Domain transformation: x [-1,1] ↦ t [0,1]
84
- t (x) = FP (1 // 2 ) * x + FP (1 // 2 )
83
+ t (x) = FP (1 // 2 ) * x + FP (1 // 2 )
85
84
86
85
function integrand ((weights, nodes))
87
86
ts = t .(nodes)
88
87
prod (weights) * f (geometry (ts... )) * differential (geometry, ts)
89
88
end
90
89
91
- return FP (1 // (2 ^ N)) .* sum (integrand, zip (weights, nodes))
90
+ return FP (1 // (2 ^ N)) .* sum (integrand, zip (weights, nodes))
92
91
end
93
92
94
93
# HAdaptiveCubature
95
94
function _integral (
96
- f,
97
- geometry,
98
- rule:: HAdaptiveCubature ,
99
- FP:: Type{T} = Float64
100
- ) where {T<: AbstractFloat }
95
+ f,
96
+ geometry,
97
+ rule:: HAdaptiveCubature ,
98
+ FP:: Type{T} = Float64
99
+ ) where {T <: AbstractFloat }
101
100
N = Meshes. paramdim (geometry)
102
101
103
102
integrand (t) = f (geometry (t... )) * differential (geometry, t)
@@ -109,7 +108,7 @@ function _integral(
109
108
# Create a wrapper that returns only the value component in those units
110
109
uintegrand (uv) = Unitful. ustrip .(integrandunits, integrand (uv))
111
110
# Integrate only the unitless values
112
- value = HCubature. hcubature (uintegrand, zeros (FP,N), ones (FP,N); rule. kwargs... )[1 ]
111
+ value = HCubature. hcubature (uintegrand, zeros (FP, N), ones (FP, N); rule. kwargs... )[1 ]
113
112
114
113
# Reapply units
115
114
return value .* integrandunits
@@ -120,32 +119,32 @@ end
120
119
# ###############################################################################
121
120
122
121
function _integral_1d (
123
- f,
124
- geometry,
125
- rule:: GaussKronrod ,
126
- FP:: Type{T} = Float64
127
- ) where {T<: AbstractFloat }
122
+ f,
123
+ geometry,
124
+ rule:: GaussKronrod ,
125
+ FP:: Type{T} = Float64
126
+ ) where {T <: AbstractFloat }
128
127
integrand (t) = f (geometry (t)) * differential (geometry, (t))
129
128
return QuadGK. quadgk (integrand, zero (FP), one (FP); rule. kwargs... )[1 ]
130
129
end
131
130
132
131
function _integral_2d (
133
- f,
134
- geometry2d,
135
- rule:: GaussKronrod ,
136
- FP:: Type{T} = Float64
137
- ) where {T<: AbstractFloat }
138
- integrand (u,v) = f (geometry2d (u,v)) * differential (geometry2d, (u,v))
139
- ∫₁ (v) = QuadGK. quadgk (u -> integrand (u,v), zero (FP), one (FP); rule. kwargs... )[1 ]
132
+ f,
133
+ geometry2d,
134
+ rule:: GaussKronrod ,
135
+ FP:: Type{T} = Float64
136
+ ) where {T <: AbstractFloat }
137
+ integrand (u, v) = f (geometry2d (u, v)) * differential (geometry2d, (u, v))
138
+ ∫₁ (v) = QuadGK. quadgk (u -> integrand (u, v), zero (FP), one (FP); rule. kwargs... )[1 ]
140
139
return QuadGK. quadgk (v -> ∫₁ (v), zero (FP), one (FP); rule. kwargs... )[1 ]
141
140
end
142
141
143
142
# Integrating volumes with GaussKronrod not supported by default
144
143
function _integral_3d (
145
- f,
146
- geometry,
147
- rule:: GaussKronrod ,
148
- FP:: Type{T} = Float64
149
- ) where {T<: AbstractFloat }
144
+ f,
145
+ geometry,
146
+ rule:: GaussKronrod ,
147
+ FP:: Type{T} = Float64
148
+ ) where {T <: AbstractFloat }
150
149
error (" Integrating this volume type with GaussKronrod not supported." )
151
150
end
0 commit comments